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Section  1:  Introduction  and  literature  review 

Predicting  blast-vehicle  interactions  is  a  daunting  task.  Typically,  extensive  experimental 
testing  is  performed,  utilizing  obsolete  or  damaged  vehicles,  thereby  precluding  an  accurate 
indication  of  blast  response  for  vehicles  currently  deployed  in  the  field.  As  an  alternative, 
numerical  simulation  can  be  employed  to  predict  these  interactions.  The  fidelity  of  the  numerical 
simulations  can  vary  from  hydrocode  simulations  to  decoupled  Computational  Fluid  Dynamics 
(CFD)  and  Computational  Structural  Dynamics  (CSD)  simulations  to  fully-coupled  fluid 
structure  CFD/CSD  simulations. 

This  report  describes  research  and  development  activities  supported  by  contract  W56HZV- 
08-C-0126,  which  was  funded  under  the  SimBRS  program,  to  develop  a  strategy  to  perform  high- 
fidelity  simulations  of  blast-vehicle  interactions,  including  complex  phenomena  such  as 
detonation,  coupled  fluid- structure  interactions,  and  non-linear  structural  response  to  severe, 
time-dependent  pressure  loads  and  debris  impact.  This  effort  leveraged  ongoing  development 
efforts  for  the  high-fidelity  flow  solver  Loci/CHEM  to  produce  a  fully  Eulerian  numerical 
capability  to  simulate  blast  phenomena  in  fluids  and  solids.  The  resulting  code  was  christened 
Loci/BLAST.  Included  in  this  report  are  the  following:  an  assessment  of  existing  technology;  a 
discussion  of  the  implementation  of  the  Jones-Wilkins-Lee  equation  of  state  and  a  prescribed 
bum  model  within  Loci/BLAST;  results  from  simulations  of  blast  interactions  with  realistic  rigid 
geometries;  a  discussion  of  the  implementation  of  the  models  needed  to  simulate  blast-soil 
interactions  within  Loci/BLAST  and  the  associated  validation  activities;  results  from  simulations 
of  blast  interactions  with  debris;  and  preliminary  results  from  coupled  fluid-structure  interaction 
simulations  using  Loci/BLAST  and  LS-DYNA.  The  remainder  of  this  section  consists  of  a 
review  of  the  state-of-the-art  in  several  blast-related  topic  areas. 

1.1  Blast  wave  interactions  with  vehicles  and  structures 

Simulations  for  simple  and  complex  geometries 

In  designing  new  armored  vehicle  concepts,  protection  against  landmine  explosions  is  of 
utmost  importance.  The  level  of  mine  protection  and  occupant  safety  largely  defines  the  vehicle 
structure  and  configuration  of  the  crew  compartment.  Therefore,  determination  of  the  pressure 
loads  applied  to  the  vehicle  structure  is  a  necessary  first  step  in  the  design  of  protective  measures 
for  vehicles. 

Computing  the  pressure  loading  and  determining  the  structural  response  is  a  demanding  task 
requiring  the  use  of  numerical  simulations.  The  numerical  simulations  can  vary  from  hydrocode 
simulations  to  decoupled  Computational  Fluid  Dynamics  (CFD)  and  Computational  Structural 
Dynamics  (CSD)  simulations  to  fully-coupled  fluid  structure  CFD/CSD  simulations. 

Numerous  experimental  and  numerical  blast  simulations  have  been  conducted  using  a 
simplified  geometry  consisting  of  a  plate  structure  [Gupta,  1999,  Gupta,  2001,  Jacinto  et  al., 
2002],  Showichen  et  al.  [Showichen  et  al.,  2005]  simulated  the  response  of  a  plate  to  an  anti-tank 
blast  using  LS-DYNA.  The  loading  was  computed  using  CONWEP  [Hyde,  1991],  which 
provides  semi-empirical  blast  loading  functions  produced  by  a  variety  of  weapon  configurations. 
Kambouchev  et  al.  [Kambouchev  et  al.,  2006,  Kambouchev  et  al.,  2007a,  Kambouchev  et  al., 
2007b]  examined  the  influence  of  the  fluid-structure  interaction  in  the  dynamic  response  of  a 
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plate  subject  to  blast  loading.  The  understanding  developed  through  these  types  of  analyses  is 
crucial  for  modeling  more  complex  structures. 

Lohner  et  al.  [Lohner  et  al.,  2008]  described  methods  for  handling  geometrically  complex, 
body-fitted,  unstructured  grids  that  contain  embedded  surfaces  and  immersed  bodies.  These 
methods  allow  for  fluid-structure  interaction,  i.e.  CSD/CFD  coupling,  using  unstructured  flow 
solvers.  Examples  include  the  interaction  of  an  explosion  with  a  generic  ship  hull  and  a  generic 
weapon  fragmentation.  Tai  et  al.  [Tai  et  al.,  2007b]  presented  a  similar  method,  referred  to  as  the 
immersed  object  method,  for  flows  around  multiple  bodies  with  overlapping,  unstructured  grids. 
Swensen  et  al.  [Swensen  et  al.,  2006]  described  a  blast  computational  framework  for  performing 
simulations  of  soil  bound  explosion  and  their  effects  on  vehicles  and  their  human  occupants.  The 
framework  couples  the  MPMICE  CFD  code  with  DYNA3D  and  LS-DYNA  finite  element  codes. 

Simulations  for  buildings  and  urban  environments 

Simulations  to  predict  blast  loadings  on  buildings  and  urban  environments  have  also  been 
conducted  where  the  buildings  are  represented  by  simple  geometric  shapes  [Smith  et  al.,  2000, 
Xu  and  Lu,  2006].  Remennikov  [Remennikov  and  Rose,  2005]  modeled  a  blast  load  in  an  urban 
environment  using  an  uncoupled  analysis  with  buildings  modeled  as  rigid  blocks  and  plates.  A 
decoupled  CFD  code  was  used  to  determine  blast  loads  on  the  buildings.  Hayhurst  et  al. 
[Hayhurst  et  al.,  1996]  simulated  a  vapor  cloud  explosion  of  an  onshore  petro-chemical  plant 
using  a  CFD  code.  The  plant  equipment  consisted  of  a  large  number  of  vessels,  pumps,  piping, 
and  other  apparatus.  The  geometrical  entities  consisted  of  simple  primitives.  In  addition,  a 
simulated  detonation  of  an  explosive  in  a  storage  building,  using  the  hydrocode  AUTODYN-3D, 
demonstrated  explosions  for  storage  areas  with  a  common  roof,  fixed  walls,  and  a  relocatable 
sliding  wall,  with  all  walls  represented  as  plates.  Rose  et  al.  [Rose  et  al.,  2005b,  Rose  et  al., 
2005a]  performed  an  explosive  blast  simulation  of  a  city  landscape  using  an  adaptive  mesh  CFD 
code  to  predict  peak  overpressure  and  impulse  along  the  rigid  geometry.  Tang  [Tang,  2007] 
conducted  similar  simulations  using  the  same  urban  setting. 

Lu  [Lu  and  Wang,  2006]  conducted  a  simulation  to  determine  the  response  of  a  multi-story 
building  to  an  above-ground  blast  using  a  coupled  approach.  Using  the  AUTODYN-3d  code,  the 
air  around  the  building  was  modeled  using  an  Eulerian  approach  and  the  building  frame  and  soil 
were  modeled  using  a  Lagrangian  approach.  Baylot  [Baylot  and  Bevins,  2007]  performed  a 
similar  simulation  and  varied  the  standoff  distance  of  the  airblast  using  the  DYNA3D  code  for 
the  structural  response  of  the  building  and  the  CTH  shock  physics  code  to  compute  the  pressure 
boundary  conditions  to  apply  to  the  finite  element  model. 

Clutter  [Clutter  and  Stahl,  2004]  performed  blast  simulations  using  the  CEBAM  hydrocode  to 
assess  facility  vulnerability.  The  vulnerability  assessments  were  conducted  to  evaluate  the 
explosion  process  in  numerous  complex  geometric  configurations.  The  scenarios  ranged  from 
settings  which  included  only  a  few  structures  to  larger  domains  such  as  urban  environments  that 
included  cooling  towers  and  surrounding  buildings,  a  car  bomb  in  a  dense  urban  setting,  and  an 
offshore  platform. 

Simulations  for  vehicles 

In  developing  mine  protection  characteristics  for  a  vehicle,  the  design  of  the  occupant 
protection  systems  can  be  performed  only  in  the  context  of  the  complete  vehicle.  However, 
notably  few  published  studies  have  addressed  the  fluid-structure  interactions  between  blast  waves 
and  a  complete  vehicle. 
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Tai  et  al.  [Tai  et  al.,  2007a]  investigated  the  interactions  and  impact  effect  between  air-blast 
waves  and  an  armored  vehicle.  The  geometry  consisted  of  a  simplified  model  of  a  six-wheeled 
armored  tank,  including  the  muzzle,  turret,  wheels,  and  hull.  The  decoupled  inviscid  CFD 
simulations  (i.e.  assuming  a  rigid  vehicle)  were  used  to  examine  the  pressure  variation  and  shock 
wave  interaction  with  the  vehicle.  It  was  found  that  the  geometry  of  the  vehicle  is  a  key  factor  in 
determining  the  intensity  of  impact  during  the  propagation  of  the  shock  waves.  Silver  [Silver, 
2006]  used  a  commercial  CFD  code  to  predict  the  overpressure  of  a  large  caliber  gun  mounted 
upon  a  simplified,  rigid  tank  geometry.  Similarly,  Kim  [Kim  and  Han,  2006]  used  a  commercial 
CFD  code  and  a  structural  response  code  (GUNBLAST)  to  predict  the  response  of  and  damage  to 
an  aircraft  wing  and  equipment  mounted  in  the  aircraft  that  were  subjected  to  repetitive  blast 
waves  from  a  gun  mounted  on  the  aircraft. 

Honlinger  et  al.  [Honlinger  et  al.,  1996]  modeled  the  occupant  compartment  of  a  6x6  vehicle, 
including  components  above  the  vehicle  floor  such  as  the  transfer  boxes,  axle  shafts,  and  tank 
and  floor  liners.  This  simulation  decoupled  the  loading  simulation  from  the  structural  simulation. 
The  time-dependent  pressure  distribution  was  calculated  using  a  2D  Euler,  explicit,  higher-order 
finite  difference  code.  This  pressure  distribution  was  then  applied  to  the  finite  element  model, 
using  LS-DYNA  to  compute  the  structural  response. 

Lottati  et  al.  [Lottati  et  al.,  1996]  used  a  coupled  CFD/CSD  approach  to  design  a  blast 
deflector  to  improve  the  survivability  of  a  tactical  vehicle  subject  to  mine  blasts  at  different 
locations  under  the  crew  cab  and  wheels.  The  pressures  were  calculated  using  the  AUGUST-3D 
Euler  CFD  code.  The  assessment  of  the  structural  response  to  the  blast  load  was  computed  using 
the  DYNA  CSD  code.  A  simplified  geometry  model  of  the  eight-wheeled  vehicle  included  the 
cab,  body,  frame,  wheels,  transmission,  and  axles. 

Fairlie  and  Bergeron  [Fairlie  and  Bergeron,  2002]  used  a  coupled  Euler-Lagrangian  approach 
with  AUTODYN  hydrocode  to  perform  mine  blast  loading  simulations  to  assess  blast  effects  on 
a  light  armor  vehicle..  The  blast  simulations  demonstrated  both  above  ground  (air  blast)  and 
buried  charges  and  were  used  to  extract  local  velocity  measurements  from  various  parts  of  the 
vehicle.  For  these  simulations,  a  complete  model  of  the  vehicle  structure,  suspension  system,  and 
wheels  was  created  in  order  to  ensure  a  realistic  response  to  a  mine  detonation  under  one  of  the 
vehicle  wheels.  The  model  also  included  other  bulky  and  heavy  components  such  as  the  engine 
gearbox,  differentials,  leaf  springs,  axles,  shock  absorbers,  wheels,  and  tires. . 

Grujicic  et  al.  [Grujicic  et.  al.,  2007a]  performed  numerical  simulations  demonstrating  the 
effects  of  a  mine  blast  on  a  commercial  truck.  The  simulations  modeled  the  interactions  between 
the  detonation-products/soil  ejecta  resulting  from  the  explosion  of  a  shallow-buried  mine  and  a 
Ford  F800  single-unit  commercial  truck.  The  simulations  were  performed  using  the  Eulerian- 
Lagrangian  coupling  options  with  the  AUTODYN  hydrocode.  An  especially  detailed  geometric 
model  of  the  truck  was  used.  The  geometric  components  included  the  body,  cabin,  transmission, 
suspension,  gear  box,  fuel  tank,  wheels,  tires,  spokes,  bumpers,  fenders,  brakes,  engine  and 
mount,  batteries,  fuel  tank,  and  other  structural  members.  The  results  showed  that:  1)  the 
kinematic  response  of  the  vehicle  to  a  landmine  detonation  and  the  amount  of  blast  momentum 
transfer  is  sensitive  to  the  proportion  of  water  in  the  sand  into  which  the  landmine  is  buried  -  for 
saturated  sand,  the  tunneling-effect  gives  rise  to  localized  damage  while  for  dry  sand,  the  damage 
is  spread  out  over  a  larger  area  -  and  2)  the  presence  of  frequency  components  in  the  initial  blast¬ 
loading  impulse,  which  match  the  vehicle’s  natural  frequency,  plays  a  significant  role  in  the 
kinematic  response  and  damage  of  the  vehicle. 

Similar  to  the  previous  reference,  Grujicic  et  al.  [Grujicic  et  al.,  2008c]  also  conducted 
simulations  to  determine  the  survivability  of  a  1994  Chevrolet  Cl 500  commercial  pick-up  truck 
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subjected  to  the  detonation  of  landmines  buried  in  different  soil  types.  These  simulations  were 
also  performed  using  the  Eulerian-Lagrangian  coupling  options  with  the  AUTODYN  hydrocode. 
Again,  a  very  detailed  model  of  the  pick-up  truck  was  used,  including  the  cabin,  bed, 
transmission,  suspension,  gear  box,  fuel  tank,  wheels,  tires,  rims,  bumpers,  fenders,  brakes, 
engine  and  mount,  batteries,  radiator,  fuel  tank,  and  other  structural  members.  The  results 
demonstrated  that  the  soil  type  and  moisture  content  affect  both  the  extent  and  spatial  distribution 
of  the  damage  and  the  kinematic  response  of  the  vehicle  and  that  the  blast  load  contains  minimal 
low-frequency  content  and  therefore  does  not  provoke  a  whole-vehicle  resonance  response  when 
subjected  to  the  mine-detonation  loading. 

Fallet  [Fallet,  2008]  conducted  a  simulation  of  an  explosion  of  a  buried  landmine  under  the 
front  wheel  of  a  civilian  pickup  truck.  The  simulation  was  performed  using  a  coupled 
Eulerian/Lagrangian  approach  with  the  HyperWorks  RADIOSS  model.  The  geometry  included 
the  complete  vehicle,  specifically  a  detailed  model  of  the  cabin,  tires,  wheels,  suspension,  and 
underbody.  The  results  of  the  simulation  were  used  to  determine  where  to  reinforce  the  vehicle 
structure  with  armor  plates  to  limit  blast  penetration  into  the  vehicle.  It  was  also  used  to  design 
energy  absorbers  for  the  vehicle  seats  to  limit  the  acceleration  level  on  the  passenger. 

Additional  references  to  simulations  of  vehicles  subjected  to  mine  blasts  are  available 
[Williams,  2002a,  Williams,  2002b,  Grujicic  et  al.,  2008,  Grujicic  et  al.,  2008b];  however,  they 
were  not  accessible  during  the  course  of  the  review. 

1.2  Generic  blast  wave  modeling 

Blast  wave  modeling  can  be  performed  using  methods  of  varying  fidelity.  Classical  method 
use  analytical  self-similar  solutions  for  a  spherical  blast  wave  originating  from  a  point  source 
[Taylor,  1950a,  Taylor,  1950b],  This  theory  can  be  expanded  to  more  practical  scenarios  by 
including  effects  such  as  the  contributions  of  the  source  mass  [Freiwald  and  Axford,  1975].  The 
ConWep  [Hyde,  1991]  software  package  can  be  used  to  provide  semi-empirical  blast  loading 
functions  produced  by  a  variety  of  weapon  configurations.  These  blast  loading  functions  can  be 
used  to  compute  structural  responses  to  various  blast  scenarios.  However,  these  prescribed 
loading  functions  cannot  predict  effects  due  to  reflections  or  obstructions  and  thus  are  generally 
applicable  only  at  sufficiently  large  distances  from  the  blast  source.  In  order  to  model  these 
effects,  a  more  detailed  Computational  Fluid  Dynamics  (CFD)  solution  is  required. 

CFD  solutions  for  blast  wave  modeling  can  range  in  detail  from  modeling  an  initial  sphere  of 
hot  dense  gas,  to  modeling  the  detonation  and  chemical  decomposition  that  occurs  within  the 
explosive  material,  and  possibly  the  after-burning  that  occurs  when  the  explosive  gases  mix  with 
the  surrounding  air.  The  simplest  of  these  models  treats  the  explosion  as  the  spontaneous  rupture 
of  an  isothermal  sphere.  When  using  an  ideal  gas  solver,  this  is  accomplished  by  defining  a 
region  of  ideal  gas  with  a  density  equivalent  to  the  explosive  material,  and  an  energy  equivalent 
to  the  energy  released  by  the  explosive.  Such  models  do  not  perform  well  in  the  near-field  region, 
but  they  can  predict  reasonably  well  the  peak  overpressure  and  impulse  at  a  suitable  distance 
from  the  explosive  charge  [Grujicic  et  al.,  2007a], 

To  obtain  a  better  near- field  prediction,  one  can  employ  the  same  isothermal  sphere 
methodology,  except  replace  the  ideal  gas  equation  of  state  with  one  that  better  represents  the 
non-idealities  of  the  explosive  gases  (i.e.,  non-idealities  of  the  air  under  high  pressure  may  also 
be  required  to  correctly  capture  near  field  physics).  Under  high  pressures,  gases  depart  from  the 
ideal  gas  model  due  to  the  increasing  influence  of  short-range  intermolecular  forces  such  as  van 
der  Waals  forces.  For  explosive  gases,  there  are  many  candidate  equations  of  state  that  can 
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describe  the  explosive  gases.  Most  are  empirically  derived  to  match  shock  front  velocity  data 
from  explosive  tests.  The  most  common  of  these  are  the  Jones-Wilkins-Lee  (JWL)  [Lee  and 
Tarver,  1980]  and  the  Becker-Kistiakowsky-Wilson  (BKW)  equations  of  state. 

The  JWL  EoS  gives  pressure  as  a  function  of  density  and  energy,  as  defined  by  the  equation: 
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where  A,  B,  R 1 ,  and  R2  are  constants  associated  with  the  material  and  pQ  is  the  density  of  the 

unexploded  explosive  solid  at  standard  temperature.  In  some  cases,  it  is  desirable  to  have  the 
temperature  of  the  exploded  gas  (e.g.,  temperature  is  required  to  model  the  secondary  burning  of 
the  explosive  products  upon  mixing  with  air).  This  requirement  poses  a  problem  for  equations  of 
state  that  describe  pressure  as  a  function  of  density  and  energy.  However,  in  such  cases,  a 
thermodynamically  consistent  temperature  can  be  derived  by  using  the  Helmholtz  free  energy 
function,  y/  -  e-Ts,  where  the  constraint  that  entropy,  s,  cannot  be  destroyed  yields  the 
relations  [Baer  and  Nunziato,  1986a] 
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Using  these  relations,  one  finds  that  the  temperature  for  the  JWL  EoS  is  given  by  [Baer  and 
Nunziato,  1986a] 
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Alternatively,  a  temperature  dependent  form  of  the  JWL  EoS  [Lee  and  Tarver,  1980]  can  be 
employed  whereby  pressure  is  expressed  as 
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The  BKW  EoS  provides  pressure  as  a  function  of  density  and  temperature  as  well  as  providing  a 
mixing  rule.  This  is  a  mixed  parameter  EoS  that  defines  pressure  using 

P  =  pRT[ \  +  Xexp(fy)l  (1-6) 

where  R  is  the  mixture  gas  constant,  [3  is  a  tuned  constant,  and  %  is  a  mixture  coefficient  given 
by  the  mixing  rule 
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where  0,  a,  k  are  tuned  constants,  while  ni  and  k  are  coefficients  of  the  gas  species  which  are 

available  through  standard  references  [Mader,  1979].  This  approach  is  advantageous  because  the 
mixture  rule,  although  simple,  still  provides  a  temperature  of  the  mixture.  In  addition,  this 
mixture  rule  correctly  reduces  mixtures  to  thermally  perfect  as  densities  decrease.  However,  a 
disadvantage  of  the  BKW  EoS  is  that  all  of  the  species  must  utilize  it. 

Mixture  rules 

When  performing  a  blast-wave  simulation  in  an  Eulerian  context,  it  is  crucial  to  determine  the 
method  for  evaluating  the  thermophysical  properties  of  the  mixture  of  materials.  For  example,  as 
the  explosive  gases  expand  into  the  surrounding  atmosphere,  the  cells  at  the  interface  between  the 
explosive  gases  and  air  will  contain  both  air  and  explosive  gases.  Similarly,  when  modeling  a 
detonation  front  through  the  explosive  material,  an  evaluation  of  a  mixture  of  EoS  will  be 
needed.  The  implementation  of  the  mixture  rule  can  affect  interface  dynamics  as  well  as  the 
numerical  robustness  of  the  solver.  There  are  several  approaches  to  consider  for  formulating  the 
mixture  EoS  from  the  component  material  EoS,  as  outlined  below. 

The  straightforward  approach  to  treating  a  mixture  of  materials  is  to  assume  that  each 
component  of  the  mixture  is  immiscible  and  that  these  components  are  in  pressure  and  thermal 
equilibrium  [Luke  and  Cinnella,  2007],  A  standard  framework  exists  to  compose  individual 
material  EoS  functions  into  a  mixture  EoS  using  these  assumptions.  However,  in  the  case  of  short 
time-scale  events  such  as  blasts,  it  is  unlikely  that  sufficient  time  elapses  so  that  thermal 
equilibration  can  occur.  Without  a  separate  energy  equation  for  each  species,  there  is  little 
justification  for  a  specific  alternative.  Nevertheless,  many  mixture  rules  that  depart  from  thermal 
equilibration  have  been  used  in  simulations  of  explosive  blasts.  Usually,  these  approaches  can  be 
justified  because  of  their  simplicity.  Otherwise,  the  proper  accounting  for  thermal  nonequilibrium 
would  require  evolving  independent  energy  equations  for  each  material  in  addition  to  specifying 
a  mixture  rule.  One  such  approach  volume-averages  the  pressure  from  each  material  equation  of 
state,  assuming  that  each  material  occupied  the  entire  volume  and  possessed  the  total  energy 
[Colella  and  Glaz,  1985].  Other  approaches  [Clutter  and  Belk,  2002]  assume  a  mass- weighted 
energy  partition  between  species  with  pressure  equilibrium.  In  the  end,  these  different  mixing 
rules  have  three  potential  impacts  on  the  simulation  of  explosive  events:  1)  some  of  these 
simplified  mixing  rules  may  improve  the  performance  of  EoS  queries,  2)  the  effective  sound 
speed  of  the  mixture  is  controlled  by  the  mixing  rule  thus  impacting  the  behavior  of  the  interface 
region,  and  3)  as  a  result  of  thermodynamic  improprieties  in  some  of  the  ad-hoc  mixing  rules, 
temperature  is  not  properly  defined.  Therefore,  these  approaches  are  not  suitable  for  considering 
temperature-dependent  models  such  as  post-detonation  burning  of  gases  released  by  the 
detonation. 

Detonation  models 

In  the  near  field  region,  accurate  modeling  of  an  explosive  event  may  require  modeling  the 
progression  of  the  detonation  wave  through  the  explosive  material.  In  some  studies  [Fiserova, 
2006]  using  numerical  models,  peak  over-pressure  and  impulse  were  found  to  be  sensitive  to  the 
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point  of  initiation  (and  by  implication  resolving  the  detonation  wave  propagation).  Generally,  the 
importance  of  modeling  these  effects  depends  upon  the  objective  and  distance  from  the  target. 

Frequently,  modeling  the  progress  of  the  detonation  front  is  more  important  for  design  of 
explosives  than  for  damage  evaluation.  If  the  detonation  phenomena  are  important,  however,  then 
there  are  several  levels  of  detail  that  can  be  employed  to  model  the  propagation  of  the  detonation 
front.  Perhaps  the  simplest  of  these  is  to  prescribe  a  known  detonation  velocity  and  force 
decomposition  of  the  explosive  material  into  gas  products  based  on  the  distance  to  the  initiation 
point.  Generally  these  models  include  some  sort  of  smoothing  parameter  whereby  the  detonation 
front  is  distributed  over  several  cells  to  improve  numerical  robustness.  Models  such  as  these  are 
implemented  in  standard  commercial  codes  such  as  LS-DYNA  and  AUTODYN. 

More  detailed  models  are  also  available  whereby  the  decomposition  reactions  are  employed 
[Menikoff,  2006,  Saurel  and  Massoni,  1998,  Tarver  et  al.,  1997],  These  models  require  a  much 
finer  mesh  resolution  than  the  prescribed  detonation  front  approach  described  earlier,  and 
generally  they  are  much  more  computationally  expensive.  For  the  applications  of  interest  to  this 
research — modeling  the  effects  of  a  relatively  generic  blast  source — these  schemes  offer  much 
more  fidelity  than  is  required,  thus  no  further  discussion  of  detailed  detonation  front  modeling 
will  be  included. 

1.3  Modeling  blast-soil  interactions 

Soil  plays  an  important  role  in  characterizing  landmine  or  IED  explosive  loads.  The  soil  can 
both  dissipate  as  well  as  focus  explosive  energy.  In  addition,  the  ejected  soil  material  can 
comprise  a  significant  fraction  of  the  impulse  loading  on  the  vehicle.  Broadly  speaking,  soil  can 
be  classified  into  two  categories:  cohesionless  (sandy)  soils  and  cohesive  (clay-based)  soils.  Each 
type  of  soil  brings  its  own  modeling  requirements.  There  have  been  significant  investigations  into 
modeling  blast-soil  interactions.  Most  of  these  models  have  employed  general  purpose  simulation 
codes  such  as  LS-DYNA  or  AUTODYN  with  soil  modeled  as  a  solid  material  with  specific 
equation  of  state  and  strength  models  [Fiserova,  2006,  Grujicic  et  al.,  2007b,  Grujicic  et  al., 
2008d,  Wu  et  al.,  2004,  Gupta,  2001,  Hlady,  2004,  Martin  and  Link,  2003,  Fairlie  and  Bergeron, 
2002,  Wang  et  al.,  2004],  In  all  of  these  models,  the  soil  material  was  treated  as  a  plastic  solid 
material  and  the  failure  model  was  based  on  a  simple  limit  on  the  tensile  strength.  The  three 
components  of  most  of  these  models  included  an  equation  of  state,  a  strength  model,  and  a  failure 
model.  Fiserova  gives  a  detailed  account  of  many  of  the  modeling  approaches  [Fiserova,  2006], 

In  most  of  the  soil  models,  a  three  phase  mixture  of  solids,  water,  and  air  was  used  to  obtain 
the  equation  of  state  [Fiserova,  2006,  Wang  et  al.,  2004,  Grujicic  et  al.,  2006],  Typically,  these 
equations  of  state  were  implemented  by  prescribing  a  look-up  table  to  define  a  density  vs. 
pressure  relationship  as  well  as  a  sound  speed  vs.  density  relationship.  In  most  of  these  models, 
values  from  lower  pressures  were  extrapolated  to  high  pressure  regimes  where  no  experimental 
data  was  available.  For  strength  models,  sophisticated  models  were  developed  based  on 
mechanistic  analogies  of  springs  and  dampers  to  describe  the  soil  micro-structure  response 
[Wang  et  al.,  2004], 

There  are  several  common  themes  in  all  of  the  above  models.  First,  they  all  make  use  of  a 
three  phase  model  to  derive  an  equation  of  state;  however,  none  of  them  represent  the  soil  as 
having  independent  velocities  for  the  different  phases  in  a  multi-phase  model.  Thus,  the 
explosive  gases  released  by  detonation  do  not  pass  through  a  cloud  of  soil  particulates,  but 
instead  escape  (along  with  some  chunks  of  soil)  as  the  soil  is  tom  open  by  the  blast.  In  addition, 
depending  on  how  the  damage  progresses  through  the  soil,  some  of  the  soil  may  erode  through 
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the  use  of  free  nodes.  These  free  nodes  represent  the  material  removed  through  failure,  and 
although  they  are  removed  from  the  simulation,  they  can  then  continue  on  their  trajectory  and 
impact  nearby  surfaces.  However,  since  these  free  nodes  do  not  occupy  volume  or  have  drag  or 
have  a  pressure  coupling  to  the  fluid  flow,  this  numerical  approach  is  not  particularly  physical, 
and  certainly  does  not  respect  the  full  complexity  of  a  more  complete  multi-phase  model. 
Because  the  models  reviewed  did  not  consider  these  effects,  we  considered  additional  literature 
regarding  models  used  in  similar  circumstances  where  dense  multi-phase  flows  experienced 
highly  energetic  events  such  as  shocks  and  explosions.  An  overview  of  this  review  is  provided  in 
the  following  section. 

Modeling  energetic  dense  multi-phase  flows 

Soil  presents  a  challenge  for  accurate  modeling  of  the  explosive  interactions  because  it  is  a 
multi-phase  mixture  composed  of  air,  water,  and  solids.  In  an  event  such  as  an  explosion,  a 
portion  of  the  soil,  in  particular  the  overburden,  becomes  a  dense,  highly  energetic,  multi-phase 
cloud  of  material  that  ejects  out  of  the  developing  crater.  We  are  interested  in  determining  what 
physical  phenomena  may  be  significant  to  modeling  such  a  dense  multi-phase  cloud  of  soil 
material.  In  this  section  we  outline  models  that  have  successfully  been  employed  in  multi-phase 
shock  interaction  problems.  While  some  of  the  soil  models  did  build  upon  a  multi-phase  model  to 
construct  the  equation  of  state,  they  did  not  consider  non-equilibrium  effects  such  as  independent 
velocities  of  the  solid  and  gas  phases,  or  the  effects  of  energy  transfer  through  particle-particle 
collisions.  However,  there  have  been  numerous  studies  on  the  solid  particles  dispersion  under 
shock  or  explosion.  We  outline  a  sample  of  these  models  in  order  to  better  understand  the  physics 
of  a  dense  cloud  of  soil  particulates  that  forms  when  a  buried  explosive  is  detonated. 

In  the  explosive  dispersal  process  of  a  packed  bed  of  solid  particles  saturated  with  explosive, 
the  flow  topology  ranges  from  a  granular  flow  during  the  propagation  of  the  detonation  within 
the  charge  to  a  dilute  gas-solid  flow  at  distances  far  from  the  source.  In  the  dilute  flow,  where 
distance  between  particles  is  great  enough,  no  pressure  terms  and  shear  stress  would  be  required 
to  describe  the  particulate  phase.  In  the  early  stage  of  the  explosion,  the  soil  particles  are  densely 
packed  and  particle-particle  collisions  have  to  be  taken  into  account. 

In  the  literature,  we  have  identified  several  proposed  two-phase  models  appropriate  for 
describing  a  dense  gas-solid,  or  a  granular  flow,  in  which  particle-particle  collisions  are 
important.  Baer  and  Nunziato  [Baer  and  Nunziato,  1986b]  (BN  model)  presented  a  two-phase 
mixture  theory  to  describe  the  deflagration-to-detonation  (DDT)  in  reactive  granular  material. 
This  theory,  which  was  based  on  a  continuum  assumption  and  included  the  effects  of 
compressibility  for  all  of  the  phases  and  compaction  for  the  granular  material,  described  a 
granular  explosive  in  terms  of  a  gas  phase  that  filled  the  interstitial  pores  between  chemically 
reacting  solid  grains.  The  velocity,  temperature,  and  pressure  of  the  two  phases  were  allowed  to 
be  unequal.  With  the  introduction  of  Helmholtz  free  energy  and  by  requiring  the  model  to  satisfy 
the  entropy  inequality  (second  law  of  thermodynamics),  specific  constitutive  equations 
representing  important  aspects  for  granular  materials  were  developed.  A  dissipation  inequality  for 
the  mixture  was  employed  to  formulate  coupling  source  terms,  requiring  that  at  each  point  the 
mixture  entropy  not  decrease  with  time.  In  the  BN  model,  the  pressure  in  the  solid  grains  equaled 
the  pressure  in  the  gas  phase  plus  the  pressure  due  to  contact  force  between  the  grains  (i.e., 
configuration  pressure)  which  depends  on  the  rate  of  change  of  Helmholtz  free  energy  with  solid 
volume  fraction.  Later,  Powers  et  al.  [Powers  et  al.,  1990]  noted  some  inconsistency  in  the 
treatment  of  the  solid  free  energy’s  volume  fraction  dependence  and  developed  a  more  general 
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model  that  incorporated  most  features  of  the  BN  model  plus  a  set  of  specific  constitutive 
relations.  Bdzil  et  al.  [Bdzil  et  al.,  1999]  examined  the  BN  mixture  model,  with  particular 
attention  paid  to  the  manner  in  which  its  constitutive  functions  were  formulated.  Deficiencies  and 
inconsistencies  in  the  derivation  in  the  BN  model  are  cited,  and  improvements  are  suggested.  It  is 
noted  that  the  entropy  inequality  constraints  do  not  uniquely  determine  the  phase  interaction 
terms.  The  resulting  flexibility  is  exploited  to  suggest  improved  forms  for  the  phase  interaction. 
The  BN  model  treats  each  phase  as  compressible  and  in  complete  thermodynamic 
nonequilibrium,  Therefore,  it  is  suited  for  description  of  the  various  thermal  and  mechanical 
processes  leading  to  detonation. 

In  contrast  to  the  BN  model,  many  two-phase  fluid  models  assume  pressure  equilibrium.  The 
volume  fraction  of  the  solid  phase  is  determined  by  an  algebraic  equation  rather  than  a  particle 
differential  equation.  In  a  conventional  two-fluid  model  for  a  dense  phase  flow,  additional  forces 
due  to  the  particle-particle  interaction  are  added  to  the  momentum  equation  for  particulate  phase. 
These  forces  include  the  pressure  and  shear  in  the  solid  phase.  Several  proposed  formulas  have 
modeled  the  solid  shear  stress  and  pressure  terms  [Harris  and  Amsden,  1994,  Gidaspow,  1994], 
Among  them,  Gidaspow’s  model  [Gidaspow,  1994],  based  on  the  application  of  kinetic  theory, 
has  drawn  significant  attention  in  the  community.  Because  of  similarities  between  particle- 
particle  interactions  and  molecular  interactions  in  a  gas,  the  concepts  from  kinetic  theory  were 
applied  to  develop  the  governing  equations  for  dense  flows:  particle -particle  collisions  are 
responsible  for  momentum  and  energy  transfer  in  a  dense  flow  in  the  same  way  that  molecular 
interactions  are  responsible  for  pressure  wave  propagation  and  viscosity  in  a  single  phase  fluid. 
Thus,  the  fluctuating  kinetic  energy  associated  with  particle-particle  collisions  is  defined  as 
granular  temperature  and  one  extra  governing  equation  for  granular  temperature  based  on  the 
Boltzmann  equation  is  added  to  the  equations  for  the  particulate  phase.  The  granular  temperature 
0  is  defined  as  one-third  mean  square  fluctuational  velocity  of  particle  motion  as  follows: 

0=|<C2>,  (1.8) 

where  C  is  the  fluctuational  velocity  of  the  particle  motion.  Granular  temperature  can  be 
produced  by  a  shearing  action  in  granular  flow  and  by  hydrodynamic  forces.  Dissipation  can 
occur  through  inelastic  particle-particle  and  particle-wall  collisions  and  dissipation  in  the  fluid. 
Granular  temperature  can  also  be  diffused  in  the  same  manner  as  heat.  The  interaction  of  the 
fluctuating  and  mean  motion  of  the  particles  gives  rise  to  an  effective  shear  viscosity  for  the 
particulate  phase,  and  particle-particle  collisions  also  generate  the  pressure  of  particle  phase. 
Both  shear  viscosity  and  pressure  of  particulate  phase  can  be  expressed  as  the  function  of 
granular  temperature,  which  provides  complete  closure  models  for  the  equation  system. 
Gidaspow’s  approach  has  the  advantage  of  being  able  to  clarify  many  physical  concepts  raised  by 
particle-particle  collisions  using  granular  temperature. 

There  are  some  heuristic  constitutive  models  for  solid  flow.  Zhang  et  al.  [Zhang  et  al.,  2001] 
estimated  the  pressure  and  sound  speed  of  the  particle  system  using  a  heuristic  interpolation 
method,  which  is  realized  by  applying  a  weighting  function  to  the  solid  volume  fraction  between 
the  solid  limit  and  the  dilute  flow  limit. 
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Models  for  shocked  multi-phase  flows 

The  literature  extensively  documents  models  that  account  for  the  interaction  between  solid 
particles  and  an  incident  shock  wave.  For  example,  Fan  et  al.  [Fan  et  al.,  2007]  investigated  both 
experimentally  and  numerically  the  interaction  of  a  planar  shock  wave  with  a  loose,  dusty  bulk 
layer,.  Experiments  were  conducted  in  a  shock  tube.  The  incident  shock  wave  velocity  and 
particle  diameters  were  measured  with  the  use  of  pressure  transducers  and  a  Malvern  particle 
sizer,  respectively.  The  flow  fields  of  both  the  gas  and  dense  particle  phase  induced  by  shock 
waves  were  visualized  by  means  of  shadowgraphs  and  pulsed  X-ray  radiography  with  particle 
tracers  added.  The  particles,  which  were  deposited  homogeneously  inside  the  test  chamber,  were 
ultra- fine  starch  particles  with  a  mean  diameter  of  15  \xm.  In  their  numerical  work,  the  Gidaspow 
model  [Gidaspow,  1994]  was  employed  to  provide  closure  for  the  particle-particle  interactions. 
Both  experimental  measurements  and  numerical  results  showed  that  interaction  of  the  shock  with 
the  granular  particle  layer  produced  curvature  in  the  incident  shock  in  the  downstream  direction. 
This  interaction  is  caused  by  the  large  acoustic  impedance  of  the  granular  material  compared  to 
that  of  the  gas.  In  addition,  under  the  action  of  incident  shock  wave,  the  transmitted  shock  wave, 
which  is  induced  inside  the  granular  material,  makes  an  oblique  angle.  Particles  are  most 
concentrated  at  the  oblique  interface  where  the  gas  and  particles  meet  behind  the  shock. 
Numerical  results  revealed  that  the  fluctuating  motion  of  the  particles,  which  dominates  the 
collision  pressure,  the  collision  momentum  transfer,  and  the  granular  conductivity,  happened 
mainly  in  the  region  behind  the  transmitted  shock  in  the  granular  material. 

Experimental  and  numerical  studies  of  shock  impingement  on  a  particle  bed  were  conducted 
in  a  vertical  shock  tube  [Rogue  et  al.,  1998].  The  shock- induced  motion  of  particle  beds  of 
various  thicknesses  was  investigated.  Analysis  of  shock  formation  showed  that  the  particle  bed 
acted  as  an  obstacle  to  the  gas  flow,  leading  to  a  reflected  shock  that  propagated  upstream,  a 
complex  refracted  shock  that  passed  through  the  bed,  and  a  weak  transmitted  shock  propagating 
downstream.  While  the  dilution  of  the  particle  bed  proceeded  (the  front  layer  of  particles  moved 
faster  than  the  lower  layer  of  particles,  leading  to  the  dilution  of  particle  bed),  a  rarefaction 
propagated  toward  the  upstream  and  a  compression  toward  the  downstream.  While  the 
rarefaction  waves  overtook  the  reflected  shock  and  weakened  it,  the  compression  waves  merged 
with  the  transmitted  waves  and  strengthened  it.  Due  to  the  stronger  blockage  effect  of  a  double 
layer  (compared  to  a  single  layer),  the  strength  of  the  reflected  shock  was  stronger  for  the  double 
layer.  In  contrast,  the  strength  of  transmitted  shock  through  the  particle  bed  was  weaker  for  a 
double  layer.  The  effect  of  inter-particle  collisions  has  been  illustrated  with  the  two-layer  bed, 
which  became  more  dilute  than  a  single  particle  layer  (the  front  layer  of  particles  moved  faster 
than  the  lower  layer  of  particles,  leading  to  the  dilution  of  particle  bed).  The  collisions  between 
particles  in  the  two-layer  bed  permitted  a  faster  dilution  of  the  particle  cloud  and  hence  presented 
less  of  an  obstacle  to  the  gas  flow.  Therefore,  the  gas  flows  faster  through  a  two-layer  particle 
bed  than  a  single  layer  one,  which  implies  that  a  single-layer  bed  presented  more  of  an  obstacle  to 
the  gas  flow.  The  dynamics  of  a  thicker  bed  (thicker  than  two  layers)  were  also  characterized  by 
faster  dilution;  however,  in  this  case,  the  rise  in  the  number  of  particle  layers  did  not  allow  a 
rapid  gas  flow  through  the  cloud  particles.  The  high  particle  density,  in  spite  of  an  important 
dilution  of  the  cloud,  still  exerted  a  significant  resistance  to  the  gas  flow.  In  their  numerical 
simulations,  Saurel  et  al.  [Saurel  et  al.,  1992]  investigated  one-dimensional  flows.  The  drag 
coefficient  used  comes  from  an  experimental  investigation  for  a  single  particle.  The  collisions 
between  particles  were  taken  into  account  by  the  use  of  a  particle  interaction  tensor,  which  is  a 
function  of  the  particle  volume  fraction.  In  their  numerical  model,  artificial  diffusion  was 
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introduced  to  handle  solution  discontinuities.  The  numerical  results  showed  that  the  drag  force  is 
the  main  factor  contributing  to  the  dynamics  of  the  dense  particle  cloud.  Inter-phase  heat  transfer 
scales  and  particle-particle  collisions  at  the  early  stage  of  cloud  motion  had  only  a  minor  effect. 

Large  solid  volume  fraction  effects  in  multi-phase  flows 

For  very  dense  multi-phase  flows,  the  effects  of  large  volume  fractions  of  the  solid  and  liquid 
phases  must  be  taken  into  account.  In  addition  to  the  particle-particle  collision  effects  outlined  in 
the  previous  section,  the  existence  of  the  large  volume  fraction  modifies  the  effective  speed  of 
sound  of  the  respective  phases.  The  speed  of  sound  drops,  sometimes  significantly,  at  large 
volume  fractions.  The  effects  of  the  change  of  sound  speed  can  be  accommodated  by  utilizing  an 
integrated  equation  of  state  that  incorporates  the  compressibility  of  each  phase.  Alternatively 
heuristic  models  can  be  used.  For  example,  Zhang  et  al.  employed  an  asymptotic  rule  [Zhang  et 
al.,  2001]  and  developed  a  heuristic  model  to  describe  the  dependency  of  particulate  phase  sound 
speed  on  the  solid  volume  fraction.  Another  challenge  is  characterizing  the  drag  coupling 
between  the  phases  when  the  volume  fraction  is  large.  The  drag  of  particles  in  dense  two-phase 
flows  differs  from  the  drag  in  dilute  flows:  The  influence  of  neighboring  particles  on  the  gas  flow 
renders  traditional  drag  correlations  invalid.  Also,  analytical  models  are  difficult  to  derive 
because  an  adequate  model  must  account  for  the  surface  of  every  nearby  particle.  Experimental 
studies  are  prevented  by  the  difficulties  of  measuring  the  force  and  the  local  flow  field  for  an 
individual  particle  in  a  cloud.  Most  of  the  data  for  particle  drag  have  been  inferred  from 
sedimentation  and  fluidization  bed  studies  in  which  drag  force  on  particles  balances  with  the 
pressure  drop  in  a  packed  bed.  Based  on  Ergun’s  pressure  drop  equation  [Ergun,  1952], 
Gidaspow  [Gidaspow,  1994]  provided  a  particle  drag  coefficient  formulation  that  characterized 
drag  in  a  two-phase  flow  with  large  solid  volume  fraction  (solid  volume  fraction  larger  than  80 
percent). 

Modeling  impulse  in  a  multi-phase  explosion 

The  damage  to  a  nearby  structure  from  a  homogeneous  explosive  is  primarily  due  to  the 
momentum  transferred  from  the  gas  flow.  In  the  near  field,  the  high  momentum  flux  of  the  gas- 
solid  flow  from  a  heterogeneous  explosive  can  generate  particularly  large  forces  on  adjacent 
bodies.  In  multi-phase  explosive  events,  the  momentum  of  the  blast  wave  is  augmented  by  the 
solid  fragment  impact.  Frost  et  al.  [Frost  et  al.,  2007]  carried  out  an  experiment  and  numerical 
work  to  characterize  the  particle  momentum  flux  generated  by  the  detonation  of  a  heterogeneous 
explosive  consisting  of  a  packed  bed  of  inert  particles  to  determine  the  relative  contribution  of 
the  particles  and  the  gas  blast  wave  to  the  impulse  exerted  on  the  nearby  structure.  In  their 
experimental  work,  the  momentum  flux  was  measured  by  the  bending  work  done  on  a  cantilever 
gage  or  the  momentum  imparted  to  a  free  piston.  In  their  numerical  model,  nonequilibrium 
processes  between  each  phase  in  the  explosive  were  considered  and  BN  [Baer  and  Nunziato, 
1986b]  model  to  represent  the  constitutive  laws  and  the  evolution  of  solid  volume  fraction  was 
adopted.  In  the  simulations,  only  drag  force  and  pressure  gradient  of  solid  phase  were  included  in 
the  particulate  phase  momentum  equation.  Other  forces  such  as  the  pressure  gradient  within  the 
gas,  the  Basset  force  (this  term  addresses  the  temporal  delay  in  boundary  layer  development  as 
the  relative  velocity  changes  with  time)  and  added  mass  force  (this  term  represents  the  force  to 
accelerate  the  gas  near  the  particle  surface)  are  neglected.  When  particles  are  added  to  a 
homogeneous  liquid  explosive,  the  peak  blast  wave  overpressure  and  the  positive  pressure 
impulse  (only  from  gas)  are  decreased  in  comparison  with  a  homogeneous  charge.  This  is  caused 
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by  the  fact  that  a  portion  of  energy  and  momentum  released  from  the  gas  goes  into  heating  and 
accelerating  the  particles.  The  total  impulse  applied  to  a  nearby  structure  by  a  particle-laden  flow 
is  the  sum  of  the  impulse  applied  by  the  gas  dynamics  and  that  from  the  particle  collision  with  the 
structure.  Comparing  the  impulse  from  the  heterogeneous  charge  (a  packed  bed  of  steel  beads 
saturated  with  sensitized  nitromethane)  with  that  of  homogeneous  charge  with  the  same  amount 
of  explosive,  the  heterogeneous  charge  produce  larger  impulse  on  the  nearby  structure  by  a  factor 
of  two  in  their  experimental  setup.  Since  the  addition  of  particles  reduces  the  gas-phase  impulse 
by  a  factor  of  two,  the  integrated  momentum  flux  from  the  particles  in  the  near  field  therefore 
must  increase  by  a  factor  of  about  four  over  the  gas  momentum.  They  concluded  that  the  total 
impulse  applied  to  a  nearby  body  was  dominated  by  the  impulse  due  to  particle  impacts  for  the 
heterogeneous  charge. 

Simulation  of  blast  interactions  with  debris 

Within  the  study  of  clouds  of  small  particles  (dusty  gas  flows),  it  is  difficult  to  ascertain  the 
interaction  effects  between  an  incident  shock  and  a  particle.  Hence,  in  experimental  studies  with 
a  single  or  small  number  of  small  particles,  it  is  difficult  to  arrive  at  the  dependency  between  the 
drag  coefficient  and  Reynold’s  number.  Due  to  refraction,  reflection,  and  diffraction,  the 
geometry  and  amplitude  of  the  wave  fronts  are  altered  due  to  the  interaction  of  the  shock  wave 
with  particles.  These  interaction  effects  become  significant  when  large  debris/fragments  are 
present  and  thus  numerous  experimental  investigations  have  been  performed,  particularly  with 
spherical  geometries.  Simulations  treating  these  large  debris/fragments  as  rigid  bodies  and 
tracking  their  interactions  with  forces  due  to  blast  waves  demonstrated  that  six  degree  of  freedom 
trajectory  computations  can  be  employed  in  the  validation  of  phenomenological  debris/ffagment 
models.  To  follow  is  a  brief  review  of  methods  for  simulating  and  tracking  moving  objects  in 
fluid  fields  and  a  list  of  exemplary  experimental  studies  by  which  the  current  models  may  be 
compared. 

Although  other  approaches  have  been  proposed,  the  two  most  appropriate  techniques  for 
simulating  and  tracking  moving  objects  in  fluid  fields  when  large  motion  is  present  consist  of 
embedded  domain  and  over-set  mesh  methodologies.  To  this  end,  the  volume  of  literature  using 
these  approaches  for  moving  body  simulations  in  aerodynamic  and  hydrodynamic  applications  is 
vast.  Representative  examples  of  embedded  approaches  may  be  found  in  [Lohner  et  al.,  1999, 
Murman  et  al.,  2003,  and  Pember  et  al.,  1995],  and  for  the  over-set  mesh  method  in  [Buning  et 
al.,  2004,  Chen  et  al.,  2008,  and  Lijewski  and  Suhs,  1994], 

The  over-set  mesh  method  overlaps  body-fitted  grids  for  all  objects  within  the  field.  The 
governing  equations  are  then  solved  on  each  grid,  and  inter-grid  communication  is  accomplished 
via  interpolation  at  boundaries  of  their  domains.  Similar  to  embedded  approaches  discussed 
below,  cells  that  fall  within  solid  bodies  may  be  deactivated.  Additionally,  since  the  inter-grid 
communication  is  typically  accomplished  at  outer  boundaries,  large  portions  of  the  grids  that  are 
overlapped  may  be  deactivated  as  well.  This  approach  requires  substantially  less  code 
complexity,  and  may  be  implemented  in  parallel  software  in  a  scalable  manner.  Search 
algorithms  are  required  to  search  adjacent  grids  in  order  to  reestablish  the  interpolation  stencils 
for  moving  bodies  as  the  simulation  proceeds.  The  body-fitted  grids  allow  a  precise  description 
of  the  geometry  and  no  issues  arise  with  regards  to  boundary  resolution  for  viscous  simulations. 
Special  care  must  be  taken  in  the  case  of  contact  between  objects,  with  regards  to  boundary 
condition  enforcement  as  well  as  potential  overlapping  between  regions  of  the  mesh  with 
disparate  grid  resolutions.  In  the  latter  case,  mesh  refinement  may  be  used  to  address  this 
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disparity.  Another  critical  issue  arises  with  topological  changes  of  a  given  body,  such  as  those 
associated  with  general  facture  of  a  structure.  The  original  structure  would  be  discretized  with  a 
given  mesh.  A  priori  details  of  the  fracture  are  not  readily  known.  This  concern  remains 
problematic  due  to  the  fact  that  each  of  the  individual  fragments/debris  would  require  new  grids 
to  be  generated.  Additionally,  each  fragment  would  essentially  be  in  contact  or  near  contact  with 
many  others. 

The  embedded  domain  method,  or  Cartesian  grid  method,  inserts  the  objects  onto  a  mesh  that 
discretizes  the  entire  domain.  A  search  algorithm  then  identifies  the  cells/edges  that  intersect  or 
fall  within  the  body.  The  cells  that  fall  within  the  body  are  deactivated.  For  those  that  intersect 
the  body,  boundary  coefficients  are  formed  to  ensure  flux  balance.  This  embedded  boundary 
treatment  is  an  approximation  and  introduces  errors  into  the  solution.  Adaptive  mesh  refinement 
may  be  used  to  reduce  these  errors  by  better  representing  the  geometric  surface.  Additionally, 
viscous  simulations,  which  require  highly  stretched  meshes  near  the  surface  for  boundary  layer 
resolution,  are  problematic.  This  issue  stems  from  the  fact  that,  since  the  mesh  is  not  body-fitted, 
regions  of  the  mesh  would  essentially  be  refined  in  nearly  an  isotropic  manner  and,  hence,  may 
result  in  prohibitively  large  grid  sizes.  Furthermore,  as  simulations  for  moving  body  applications 
proceed,  the  geometric  description  would  need  to  be  retained,  while  the  search  algorithm  used  to 
re-identify  the  cells  that  intersect  and  fall  within  the  body,  extrapolation  of  the  solution  for  cells 
that  re-enter  the  domain  (i.e.,  deactivated  points  become  part  of  the  computational  domain),  and 
adaptive  refinement  all  performed  again.  However,  topological  changes  and  general  fracture  of  a 
structure  does  not  pose  any  difficulties  within  the  embedded  domain  method. 

There  are  numerous  experimental  investigations  concerned  with  studying  the  interaction 
between  large  particles  and  shock  waves.  A  short  list  of  potential  studies  for  comparison  includes 
[Britan  et  al.,  1995,  Ingra  and  Takayama,  1993,  and  Susuki  et  al.,  2005],  These  studies  consisted 
of  utilizing  shock  tubes  to  accelerate  spherical  particles  of  various  diameters  and  densities. 
Within  these  experimental  investigations,  the  trajectories  and  velocities  of  the  particles  were 
measured.  From  these,  calculations  of  the  drag  coefficient  were  made  and  typically  presented  as  a 
function  of  the  relative  Reynold’s  number.  Additionally,  in  Susuki  et  al.  2005,  due  to  large 
vertical  components  of  the  trajectories,  the  rotational  speed  of  the  particles  was  also  measured  to 
ascertain  the  influence  of  the  Magnus  effect.  Although  large  rotational  speeds  were  measured,  it 
was  concluded  that  the  Magnus  force  had  little  influence  and  did  not  significantly  contribute  to 
the  lift  of  the  particles.  Studies  performed  in  smooth  and  rough  wall  horizontal  shock  tubes 
indicated  that  the  floor  conditions  had  the  strongest  influence. 

1.4  Blast  generated  fluid-structure  interaction  (FSI) 

Generally,  there  are  two  approaches  by  which  the  structure  is  encompassed  within  the 
surrounding  fluid  media.  One  approach  represents  a  non-overlapping  technique  whereby  the  CFD 
and  CSD  domains  share  a  common  wetted  surface  (outer  mold  line).  In  this  approach  the  surface 
discretizations  and  corresponding  volume  meshes  must  be  adapted  to  reflect  the  structural 
deformations.  MSU  currently  has  well-validated  software  to  accommodate  this  methodology 
[Blades  and  Newman,  III,  2007a, b]  for  FSI  of  aerospace  vehicles.  It  is  recognized  that  this 
approach  is  computationally  the  most  efficient,  but  it  has  potential  bottle  necks  for  large 
structural  deformations  that  could  cause  a  failure  in  the  volume  adaptation  process  as  well  as  an 
inability  to  accommodate  geometry  fragmentation.  The  other  approach  is  an  overlapping 
(embedded)  technique  whereby  the  CSD  mesh  is  overlapped  onto  the  CFD  domain.  Although 
more  algorithmically  complex  and  computationally  expensive,  this  approach  suffers  no  such 
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shortcomings  with  regard  to  geometry  fragmentation.  However,  the  interdisciplinary  transfer  of 
data  in  the  cases  of  surface  cracking  or  fragmentation  may  require  solution  adaptive  meshing 
techniques  in  order  to  ensure  adequate  resolution.  Technology  for  this  approach  is  currently 
under  development  within  the  Loci/CHEM  framework/code. 

Newman  [Newman,  III,  1997]  provided  an  early,  comprehensive  review  of  non-overlapping 
approaches  for  FSI;  his  effort  will  not  be  duplicated  here.  More  recently,  Jaiman  et  al.  [Jaiman,  et 
al.,  2006]  further  classified  these  non-overlapping  approaches  into  point-to-point,  point-to- 
element,  and  common-refinement  schemes.  Within  this  context,  the  current  MSU  technology 
falls  into  the  point-to-point  category.  Jaiman  then  additionally  subdivided  the  point-to-element 
schemes  into  node-projection  [Farhat,  et  al.,  1998]  and  quadrature-projection  [Cebral  and  Lohner, 
1997]  methods.  Jaiman  compared  the  two  traditional  point-to-element  methods  with  their 
developed  common-refinement  scheme  for  a  variety  of  FSI  applications.  Both  point-to-element 
approaches  illustrated  an  oscillatory  solution  for  interface  displacements  and  shear  stress  for  blast 
wave-solid  interaction.  Furthermore,  errors  in  excess  of  20%  were  present  when  compared  to 
results  from  matching  meshes.  Within  the  study  of  [Jaiman,  et  al.,  2006],  point-to-point  schemes 
were  not  considered  and,  hence,  will  require  investigation  for  their  applicability  and  accuracy  to 
blast  wave-structural  interactions. 

Although  an  exhaustive  list  will  not  be  given  here,  there  are  many  researchers  currently  using 
loosely  coupled  CFD/CSD  methodologies  to  simulate  blast-structure  interaction  problems. 
Probably  the  most  notable  contributions  have  been  given  by  Baum,  Lohner,  and  their  co-workers 
[Baum,  et  al.,  2003,  Lohner,  et  al.,  2004b, c,  Baum,  et  al.,  2004a,  Lohner,  2004a,  Baum,  et  al., 
2006,  Rice,  et  al.,  2006,  Baum  and  Lohner,  2006,  Soto,  et  al.,  2008,  Baum,  et  al.,  2008],  In  these 
references,  the  Euler/Reynold’s  Averaged  Navier-Stokes  solver  FEFL098  [Lohner  and  Baum, 
1992]  and  MARS3D  [Pelessone  and  Charman,  1998]  (a  derivative  of  DYNA3D  [Whirley  and 
Hallquist,  1991])  were  utilized  for  the  CFD  and  CSD  analyses,  respectively.  Early  work  by  these 
researchers  utilized  the  non-overlapping  quadrature-projection  scheme  of  [Cebral  and  Lohner, 
1997],  but  later  adopted  the  embedded  approach  [Baum,  et  al.,  2004b],  A  variety  of  blast  wave- 
structure  interactions  were  analyzed  within  this  body  of  work.  These  simulations  ranged  from 
blast  loading  on  plates  to  building  structures.  Furthermore,  the  fidelity  of  the  simulations  also 
ranged  from  relatively  coarse  grid  Euler  simulations  to  detailed  analyses  of  blast  wave  evolution. 
Of  particular  interest  to  the  current  work  are  Baum,  et  al.,  2006  and  Baum,  et  al.,  2008,  which 
dealt  with  the  blast-wave  structural  interaction  on  plates.  For  example,  in  Baum,  et  al.,  2006  the 
numerical  solutions  of  bare  charges  within  a  cylindrical  tunnel  were  directly  compared  with 
experimental  results.  Pressure-time  histories  were  compared  for  blast  wave  impingement  on  rigid 
end  walls  and  deformations  compared  for  deformable  end  plates.  Additional  examples  of  the 
loosely  coupled  CFD/CSD  methodology,  as  applied  to  plates  and  shells,  may  be  found  in  [Liang 
and  Hsu,  2001]  and  [Chan,  2004],  In  [Liang  and  Hsu,  2001]  unsteady  cylindrical  blast- wave 
interaction  with  a  flat  plate  was  simulated,  whereas  in  [Chan,  2004]  the  FSI  of  a  flexible  tent 
structure  was  modeled. 

As  directed  by  the  COR,  the  FSI  coupling  shall  utilize  the  LS-DYNA  code  [LS-DYNA  User 
Manual,  2007]  for  the  CSD  analyses  in  this  study.  LS-DYNA  is  a  commercial  code,  and  currently 
no  API  exists  for  coupling  with  general  purpose  CFD  codes.  Hence,  coupling  must  be 
accomplished  via  the  development  of  an  I/O  file  management  system.  Similar-type  management 
systems  have  been  developed  in  the  past  for  this  purpose.  For  example,  NASA’s  FIDO 
framework  [Weston,  et  al.,  1994],  as  well  as  a  system  developed  at  Boeing  [ Borland,  1990],  was 
used  to  couple  govemment/industrial  CFD  codes  with  commercial  CSD  software  such  as 
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NASTRAN.  More  recently,  CFDRC  developed  the  commercial  package  MDICE  [Sheta,  et  al., 
1999]  for  managing  data  produced  from  multidisciplinary  analyses  involving  various  single¬ 
discipline  software  packages.  Within  the  scope  of  this  task,  although  representing  an  added 
computational  expense,  this  approach  is  feasible  for  the  demonstrations  to  assess  the  FSI 
capabilities  to  be  delivered.  Additionally,  for  one-way  coupling,  this  procedure  should  be  very 
useful. 

The  developed  coupled  air-blast-wave-structure  interaction  analysis  capability  will  be 
evaluated  with  simplified  components  and  numerical  results  compared  with  experimental  data. 
These  simulations  shall  include  analyzing  the  shock-structure  interaction  of  bare  explosives  of 
different  charge  weights  at  various  stand-off  distances  as  it  impinges  on  plates.  The  utilization  of 
plates  for  this  purpose  is  three  fold.  First,  it  provides  a  reduction  in  geometric  complexity.  This 
reduction  facilitates  the  use  of  the  non-overlapping  approach  to  FSI  (in  the  absence  of 
fragmentation),  as  well  as  lowering  the  computational  requirements  associated  with  data 
management  between  the  single-discipline  software  components.  Second,  vast  amounts  of 
experimental  data  (briefly  reviewed  below)  have  been  produced  for  the  purpose  of  preliminary 
validation  of  blast-wave  FSI.  Third,  plates  are  one  of  the  most  basic  elements  in  structures  and 
thus  are  representative  of  the  FSI  that  will  be  simulated  in  future  blast  wave-vehicle  interaction 
analyses. 

Numerous  experimental  results  are  present  in  the  literature  for  metallic,  composite,  and 
cellular  solid  plates.  The  current  review  is  limited  to  experimental  data  obtained  for  metallic 
plates  only.  Within  this  body  of  literature,  recorded  data  takes  the  form  of  measured  pressures  on 
both  the  plate  and  the  support  structures,  surface  accelerations  (particularly  when  used  to  validate 
transient  dynamic  analyses),  and  structural  deflections  during  and  after  the  blast  event.  Jacinto  et 
al.,  [Jacinto,  et  al.,  2001]  presented  an  experimental  and  computational  study  of  non-stiffened 
metallic  (steel)  plates  subjected  to  air  blast  loadings.  Within  this  data  set,  four  different  tests  with 
various  amounts  of  explosive  were  carried  out  on  two  plates.  One  plate  was  clamped  at  the  base 
while  the  other  was  clamped  at  the  four  edges.  Experimental  pressure-time  and  acceleration-time 
histories  were  recorded.  Trzcinski  and  Cudzo  [Trzcinski  and  Cudzo,  2004]  experimentally 
studied  edge-clamped  rectangular  steel  and  steel-composite  plates  subjected  to  air  blast  waves 
generated  by  cylindrical  TNT  charges  of  various  masses.  Available  data  for  this  series  of 
experiments  is  comprised  of  resultant  force-time,  displacement-time  at  numerous  plate  locations, 
and  final  displacment  distributions  along  the  plate  center  line.  Bonorchis  and  Nurick  [Bonorchis 
and  Nurick,  2009]  ran  a  series  of  blast  loading  experiments  on  rectangular  plates  to  ascertain  the 
influence  on  the  experimental  set-up  and  demonstrate  its  importance  on  subsequent  numerical 
simulations.  Additionally,  Langdon  and  Suhleyer  [Langdon  and  Suhleyer,  2005]  conducted 
experimental  investigations  aimed  at  studying  the  influence  of  connection  details  whereby  large 
plastic  deformations  were  produced  in  the  panels  without  rupture.  Similar  experimental  and 
computational  investigations  were  presented  in  [Louca  and  Pan,  1998]  and  [Pan  and  Louca, 
1999],  In  these  articles,  stiffened  and  unstiffened  plates  were  studied,  and  the  data  presented 
included  displacement-time,  acceleration-time,  and  strain-time  histories.  Neuberger  et  al., 
[Neuberger,  et  al.,  2009]  conducted  a  series  of  experiments  in  which  the  peak  transient  dynamic 
and  the  final  residual  deflection  were  measured  in  order  to  study  the  elastic  springback  of 
metallic  plates.  Moreover,  Galiev  [Galiev,  1996]  presented  counter-intuitive  behavior  of  plates. 
Experiments  were  conducted  that  demonstrated  the  appearance  of  final  deflections  that  were 
contrary  to  the  direction  of  the  impulsive  loading.  An  explanation  of  this  behavior  was  developed 
using  theoretical  techniques  and  aerodynamic/structural  dynamic  models. 
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Section  2:  Simulation  of  blast  interactions  with 
realistic  rigid  geometries 

2.1  Geometric  model  definition 

The  objective  of  this  effor  is  to  perform  blast  simulations  on  a  complex,  realistic  vehicle 
using  the  Loci/CHEM  code  where  the  vehicle  is  assumed  to  be  rigid.  To  perform  said  blast 
simulations  required  a  sufficiently  detailed  geometric  definition  of  the  vehicle,  including 
underbody  components,  tires,  wheels,  handles,  and  mirrors,  that  also  adequately  represented  a 
military  High  Mobility  Multi-purpose  Wheeled  Vehicle  (HMMWV  or  Humvee).  Mississippi 
State  University  (MSU)  does  not  possess  a  geometric  definition  of  a  HMMWV  type  vehicle,  so 
multiple  web  sites  were  searched  to  find  potential  candidate  models.  Neither  computational  fluid 
dynamics  (CFD)  models  or  finite  element  models  (FEM)  suitable  for  analysis  were  found.  The 
initial  model  was  obtained  from  [3dcadbrowser.com_a],  However,  this  model  was  deemed  to 
have  insufficient  underbody  detail.  A  second  HMMWV  model  with  much  more  geometric  detail 
was  obtained  from  [the3dstudio.com].  This  model,  however,  was  not  approved  by  TARDEC. 

It  was  then  decided  that  a  commercial  sport  utility  vehicle  (SUV)  would  be  used  for  the  blast 
simulation.  A  Nissan  XTrail  4x4  SUV  geometry  model  was  obtained  from 
[3dcadbrowser.com_b],  which  TARDEC  approved  for  use  on  December  17,  2008.  The  SUV 
geometric  model  definition  is  discrete  (i.e.  composed  of  polygons  and  vertices),  rather  than  a 
mathematical  model  like  a  NURBS  definition.  The  discrete  SUV  geometric  model  consists  of 
204,062  triangles  and  108,372  vertices  (Figure  2.1).  The  underbody  definition  is  not  complete. 
There  are  no  surfaces  for  the  wheel  wells  or  the  section  underneath  the  engine  compartment.  In 
addition,  the  suspension  and  exhaust  components  are  floating  in  space,  dettached  from  the  rest  of 
the  vehicle. 

2.2  Geometric  model  creation 

The  original  approach  to  generate  the  volume  mesh  used  the  discrete  model  as  the  geometry 
definition  to  build  the  CFD  surface  mesh.  SolidMesh,  MSU’s  in-house-developed  CAD  and 
unstructured  mesh  generation  software,  has  the  capability  to  use  a  discrete  model  as  the 
underlying  geometric  definition.  The  discrete  geometry  definition  can  then  be  re-meshed  with  an 
appropriate  resolution  suitable  for  a  CFD  mesh.  In  order  for  the  discrete  geometry  capability  to 
work  correctly,  the  discrete  model  must  be  topologically  valid  or  watertight.  That  is,  all  of  the 
components  must  connect  to  one  another.  However,  the  chosen  discrete  SUV  model  is  not 
topologically  valid  and  thus  is  not  suitable  for  use  as  a  discrete  geometry  definition.  As  illustrated 
in  Figure  2.2,  the  triangles  representing  the  window  trim  do  not  match  those  representing  the 
window  glass  or  door  panel,  leaving  the  window  trim  floating  in  space.  In  fact,  none  of  the 
models  obtained  from  either  the  3DCADBrowser  or  3DStudio  web  sites  were  topologically 
valid.This  situation  appears  to  be  typical  for  models  of  this  type  available  on  the  web.  Their 
purpose  is  for  visualization,  which  does  not  require  models  that  are  topologically  valid,  rather 
than  for  analysis.. 


16 


UNCLASSIFIED 


(b)  Underside  view 

Figure  2.1:  Discrete  geometry  definition  of  the  Nissan  XTrail  sport  utility  vehicle. 
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Figure  2.2:  Topology  problems  with  the  discrete  model  definition. 

The  grid  repair  toolkit  within  SolidMesh  was  unable  to  repair  all  of  the  problems  with  the 
model,  requiring  significant  manual  labor  to  rectify  the  indentified  issues.  Therefore,  the  discrete 
SUV  model  was  treated  as  a  cloud  of  points.  Although  labor  and  time-intensive  as  well, 
geometric  curves  and  surfaces  were  then  manually  created  from  the  cloud  of  points  using  curve 
and  surface  generation  CAD  tools  within  SolidMesh.  NURBS  curves  were  generated  using 
splines,  and  NURBS  surfaces  were  generated  from  the  curves  using  multiple  surface  generation 
techniques,  including  transfinite  interpolation,  lofting,  and  ruled  surfaces.  A  half-plane  geometric 
model  of  the  Nissan  XTrail  4x4  sport  utility  vehicle  (SUV)  was  constructed  using  this  approach. 
Every  effort  was  made  to  retain  as  much  geometric  detail  as  possible.  A  comparison  of  the  points 
from  the  original  discrete  model  and  the  manually  generated  NURBS  surfaces  is  shown  in  Figure 
2.3.  The  white  points  shown  in  these  images  represent  the  vertices  from  the  discrete  model.  As 
illustrated,  the  geometric  CAD  surfaces  pass  through  the  points  of  the  discrete  representation  and 
thus  represent  the  surface  accurately. 

The  completed  model  is  shown  in  Figure  2.4(a)  and  details  of  the  front  and  rear  underbody 
components  are  shown  in  Figures  2.4(b)  and  2.4(c).  Detailed  views  of  the  wheel  and  suspension 
are  shown  in  Figures  2.5  and  2.6,  respectively.  The  underbody  definition  for  the  original  discrete 
model  was  not  complete.  There  was  no  surface  definition  for  the  wheel  wells  and  the  section 
underneath  the  engine  compartment.  In  lieu  of  these  missing  surfaces,  reasonable  approximations 
were  made  for  these  components.  In  order  to  create  a  topologically-valid,  watertight  model,  a 
ground  plane  and  hemispherical  outer  boundary  were  added  to  the  domain.  The  diameter  of  the 
outer  boundary  surface  was  40  meters,  which  is  approximately  10  body  lengths. 
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(a)  Front  view 


(b)  Rear  view 

Figure  2.3:  Comparison  of  the  CAD  surfaces  with  the  original  discrete  model  definition. 
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(a)  Entire  vehicle 


(b)  Front  underbody  components  (c)  Rear  underbody  components 

Figure  2.4:  Geometric  model  of  the  Nissan  XT rail  sport  utility  vehicle. 
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ciam. 


Figure  2.5:  Geometric  model  of  the  tire,  wheel,  and  rim. 


Figure  2.6:  Geometric  model  of  the  suspension. 
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2.3  Mesh  generation 

After  the  completion  of  the  topologically  valid  geometric  model,  the  next  step  involved 
suitably  discretizing  the  model  for  the  blast  simulation.  To  discretize  the  model  in  SolidMesh, 
grid  point  spacings  were  applied  to  all  of  the  geometric  vertices.  The  point  spacing  values  were 
used  to  control  the  local  element  size.  The  nominal  point  spacing  on  the  vehicle  body  wass 
1.25cm.  Some  regions  required  a  locally  smaller  point  spacing  to  adequately  represent  the 
geometric  detail.  For  the  suspension,  wheels,  brush  guard,  underbody  and  other  geometric  details, 
a  nominal  spacing  of  3.125mm  was  used.  The  global  surface  grid  was  then  generated  surface  by 
surface. 

Mesh  generation  using  SolidMesh  is  an  automated  process  using  advancing  front/local 
reconnection  algorithms.  The  edge  grid  generation  is  automatically  generated  using  a  one¬ 
dimensional  advancing  front  along  the  curves  advancing  from  the  geometric  vertices.  The  surface 
mesh  generation  is  automatically  generated  using  a  two-dimensional  advancing  front  advancing 
from  the  surface  boundaries  (i.e.,  the  curves).  The  volume  mesh  generation  is  automatically 
generated  using  a  three-dimensional  advancing  front  advancing  from  the  surfaces. 

The  geometric  model  shown  in  Figure  2.4  consists  of  1,834  geometric  surfaces.  In  order  to 
create  a  watertight  or  topologically  valid  geometric  model,  some  of  the  created  surfaces  had 
highly  acute  angles,  i.e.,  sliver  surfaces,  which  can  negatively  impact  the  surface  mesh  quality 
and  present  difficulties  for  the  flow  solver.  Therefore,  the  composite  surface  capability  within 
SolidMesh  was  used  to  generate  a  surface  mesh  of  suitable  quality.  A  composite  surface  is  a 
collection  of  surfaces  that,  for  mesh  generation  purposes,  is  treated  as  a  single  topological  entity. 

An  example  of  a  sliver  surface  is  illustrated  in  Figure  2.7(a).  This  highlighted  region,  isolated 
and  magnified  in  Figure  2.7(b),  is  part  of  the  vehicle  trim  molding  on  the  front  quarter  panel.  It 
consists  of  nine  surfaces,  for  which  the  mesh  is  shown  in  Figure  2.7(c).  As  illustrated,  the  mesh 
quality  is  poor  due  to  the  surface  topology.  A  composite  surface  was  created  using  these  nine 
component  surfaces,  and  the  resulting  surface  mesh  is  shown  in  Figure  2.7(d).  Without  the 
constraints  of  the  interior  geometric  edges,  the  composite  surface  mesh  is  more  isotropic  and 
results  in  a  better  quality  mesh.  The  minimum  angle  in  the  composite  surface  mesh  is  increased 
from  7.2  degrees  to  30.7  degrees  and  the  maximum  angle  decreased  from  122.5  degrees  to  98.1 
degrees.  Significantly,  there  is  no  approximation  in  the  final  composite  surface  mesh,  i.e.,  the 
surface  grid  points  for  the  composite  surface  lie  on  the  original  geometric  surface.  The  composite 
surface  does  not  modify  the  underlying  geometric  definitions  of  a  geometric  surface,  unlike  a 
carpet  surface  which  projects  and  interpolates  the  data.  A  total  of  76  composite  surfaces, 
consisting  of  383  component  surfaces,  were  defined.  The  locations  of  the  composite  surfaces 
(highlighted  in  yellow)  are  shown  in  Figure  2.8. 

For  the  blast  simulation,  an  isotropic  (i.e.  no  boundary-layer  resolution)  unstructured  volume 
mesh  is  suitable.  In  order  to  decrease  the  volume  grid  generation  time,  a  volume  mesh  was 
generated  for  the  half-plane  geometric  model.  This  half-plane  volume  mesh  was  then  reflected 
about  the  vehicle  plane  of  symmetry  to  obtain  the  volume  mesh  of  the  full  vehicle.  The  final  SUV 
surface  mesh  is  shown  in  figure  2.9  and  consists  of  626,843  nodes  and  1,253,858  triangles.  The 
full  vehicle  volume  mesh  contains  4,765,060  nodes  and  26,429,880  tetrahedra.  The  minimum 
dihedral  angle  of  the  volume  mesh  =  9.75  degrees  and  the  maximum  dihedral  angle  =  159.7 
degrees.  To  provide  a  sample  of  the  field  resolution,  cutting  planes  were  taken  at  various 
locations  throughout  the  volume  mesh  (Figure  2.10).  A  cutting  plane  taken  along  the  vehicle 
symmetry  plane  is  shown  in  figure  2.10(a),  and  cutting  planes  taken  near  the  suspension,  brush 
guard,  and  side  view  mirror  are  shown  in  figures  2.10(b),  2.10(c),  and  2.10(d),  respectively. 
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(a)  Location  of  component  surfaces 


(b)  Detail  of  component  surfaces 


(c)  Original  surface  mesh  (d)  Composite  surface  mesh 

Figure  2.7:  Composite  surface  example  for  surfaces  with  poor  topology. 
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Figure  2.8:  Composite  surface  definitions  for  the  Nissan  XTrail  sport  utility  vehicle. 


Figure  2.9:  Surface  mesh  for  the  Nissan  XTrail  sport  utility  vehicle. 
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(a)  Cutting  plane  along  vehicle  centerline 


(b)  Cutting  plane  taken  near  the  front  suspension 
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(c)  Cutting  plane  near  the  vehicle  brush  guard 


(d)  Cutting  plane  near  the  side  view  mirror 


Figure  2.10:  Volume  grid  cutting  planes. 
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2.4  Blast  simulations  for  rigid  vehicle 

Our  previous  validation  results  indicate  that  0(1)  mm  mesh  spacings  are  necessary  to  capture 
the  initial  phases  of  the  detonation.  However,  this  spacing  can  be  relaxed  significantly  as  the  blast 
wave  evolves.  To  reduce  computational  costs,  the  approach  we  employed  for  the  blast-vehicle 
interaction  initiated  the  blast  simulation  for  a  cylindrical  charge  as  an  axisymmetric  computation 
and  ran  this  simulation  until  the  outer  shock  wave  had  propagated  a  suitable  distance.  At  this 
point,  the  axisymmetric  blast  solution  was  interpolated  onto  the  SUV  mesh  and  the  computation 
initiated. 

WD0019  does  not  specify  the  number  of  cases  to  be  simulated  nor  does  it  specify  the  position, 
mass,  or  composition  of  the  charge.  In  an  email  (dated  June  18,  2009),  TARDEC  personnel 
suggested  that  we  use  15kg  C-4.  The  C-4  charge  was  defined  to  be  a  cylinder  with  a  diameter  of 
35.76418cm  and  a  height  of  8.9410cm.  Two  cases  were  simulated  using  different  charge 
locations.  In  the  first  case,  hereafter  refered  to  as  the  “undertire”  case,  the  geometric  center  of  the 
bottom  face  of  the  charge  was  located  on  the  ground  plane  31.558cm  ahead  of  the  front  axle  and 
2.8379cm  to  the  right  of  the  tire  symmetry  plane  at  the  point  labeled  “A,”  as  shown  in  Figure 
2.11.  The  second  case,  denoted  as  the  “centerline”  case,  had  the  charge  placed  under  the 
occupants  on  the  vehicle  centerline  appoximately  98.525cm  behind  the  axle.  In  the  “centerline” 
case,  the  vehicle  plane  of  symmetry  was  used  as  a  plane  of  symmetry  in  the  computational 
domain  in  order  to  reduce  the  cost.  Each  simulation  was  run  until  a  physical  time  of  lOOOps 
elapsed.  In  each  case,  detonation  was  initiated  at  the  center  of  the  charge.  The  HLLE  flux 
definition  was  used  [Einfeldt,  1988].  An  explicit,  second-order  Runge-Kutta  time  integration 
algorithm  was  employed  with  the  time  step  defined  based  on  a  maximum  permitted  CFL  number 
of  0.75.  The  approach  employed  to  simulate  the  detonation  and  subsequent  blast  is  discussed  in 
detail  in  Section  3. 


Figure  2.11:  Position  of  15kg  C-4  charge  for  “undertire”  case. 
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Two  sets  of  animations  were  generated,  one  for  each  charge  position.  Many  of  the  phenomena 
visible  in  the  animations  are  also  visible  in  static  images.  Figures  2.12-2.14  show  sequences  of 
images  generated  from  the  “undertire”  simulation  at  times  t=207ps,  407ps,  607ps,  807ps,  and 
1007ps.  Figure  2.12  illustrates  the  evolution  of  an  isosurface  of  the  JWL  gas  mass  fraction 
colored  by  local  pressure  (in  units  of  Pa).  This  variable  was  chosen  for  the  isosurface  because  it 
provides  a  good  approximation  to  the  blast  front.  The  impact  of  the  cylindrical  shape  of  the 
charge  on  the  blast  front  is  evident  in  Figure  2.12(a).  The  complexity  of  the  flow  interactions 
with  the  wheels  and  wheel  wells  is  apparent  in  Figures  2.12(b)-2.12(e).  Figure  2.13  shows  the 
vehicle  and  the  ground  plane  shaded  by  local  pressure.  Numerous  complex  wave  patterns  due  to 
wave  interactions  and  flow  interations  with  the  geometry  are  present  on  the  vehicle  and  the 
ground  plane.  Of  interest  are  the  diffraction  patterns  shown  in  Figures  2.13(c)  and  2.13(d)  due  to 
the  interaction  of  shocks  with  the  front  tires.  Figure  2.14  shows  the  vehicle  underbody  shaded  by 
local  pressure.  Again,  the  numerous  wave  interactions  illustrate  the  complexity  of  the  flow 
underneath  the  vehicle.  This  is  particularly  evident  in  figure  2.14(c),  which  clearly  shows  several 
shock-shock  interactions.  By  approximately  1  OOOps  from  the  initiation  of  the  blast,  the  pressure 
over  most  of  underbody  shown  in  the  image  is  approaching  the  original  pressure. 

Figures  2.15-2.17  show  a  similar  set  of  plots  for  the  “centerline”  case.  The  complexity  of  the 
resulting  flow  is  evident  from  the  figures.  Figure  2.16  shows  complex  diffraction  patterns  in  the 
ground  plane  due  to  the  interaction  of  shock  waves  with  the  tires.  Complex  wave  interactions  are 
again  present  on  the  underbody,  as  shown  in  figure  2.17. 
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(e)  t=1007  jus 


Figure  2.12:  “Undertire”  blast  simulation  -  Isosurface  of  JWL  gas  mass  fraction  colored  by 

local  pressure. 
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(a)  t=207 jus 


(c)  t=607us 


t=807jis 


(e)  t=1007|us 


Figure  2.13:  “Undertire”  blast  simulation  -  Surfaces  shaded  by  local  pressure. 
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(a)  t=207 jus 


t=807 jus 


(e)  t=1007  jus 


Figure  2.14:  “Undertire”  blast  simulation  -  Surfaces  shaded  by  local  pressure. 
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(a)  t=207us 


(b)  t=407us 


t=807jas 


(e)  t=1007|us 


Figure  2.15:  “Centerline”  blast  simulation  -  Isosurface  of  JWL  gas  mass  fraction  colored 

by  local  pressure. 
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(e)  t=1007|us 


Figure  2.16:  “Centerline”  blast  simulation  -  Surfaces  shaded  by  local  pressure. 
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(c)  t=607us 


t=807jLLS 


(e)  t=1007|us 


Figure  2.17:  “Centerline”  blast  simulation  -  Surfaces  shaded  by  local  pressure. 
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Section  3:  Implementation  of  models  to  simulate  blast- 
soil  interactions 


3.1  Governing  equations  for  the  soil/blast  modeling 

The  governing  equations  for  our  soil-blast  model  is  that  of  a  multicomponent  inviscid  flow 
assuming  thermal,  pressure,  and  velocity  equilibriation  between  all  materials  in  any  given 
location  of  space.  These  equations  are  given  by  the  conservation  of  component  mass, 
momentum,  and  energy.  The  conservation  of  component  mass  is  given  by  the  equation 

^  +  V •(p,u)=  0, Vi  s{l...NS},  (3.1) 

dt 

where  pt  is  the  component  density  of  material  i ,  and  ii  is  the  material  velocity  vector.  Note  that 
the  overall  material  density  at  any  point  in  space  is  given  by  p  =  pt .  The  momemtum 
conservation  equation  is  given  by 

—  +  V  •  {puu)  =  Vp ,  (3.2) 

dt 

where  p  is  the  pressure  of  the  component  mixture.  The  conservation  of  energy  is  given  by  the 
equation 

+  V  •  [(pe0  +  p)u]  =  0 ,  (3.3) 


where  e0  is  the  total  of  fluid  kinetic  energy  and  internal  energy  given  by  the  expression 


1  -  - 

e0  =  2W'M+eint  ernaf 


(3.4) 


These  equations  are  closed  by  a  multicomponent  equation  of  state  that  relates  pressure  to  the 
material  densities  of  the  components  and  the  internal  energy  as  represented  by 
p  -  P{pi, p2>' " > Pns ’emtemai) •  The  formulation  of  this  equation-of-state  will  be  central  to  the 


soil  and  blast  gas  models,  They  will  be  derived  from  a  mixture  rule  that  combines  single 
component  equations  of  state  for  each  material,  assuming  pressure  and  thermal  equilibriation. 


3.2  A  multicomponent  equation  of  state  for  blast  modeling 

To  support  non-ideal  equations  of  state  for  blast  modeling  using  the  explicit  time  integrator, 
we  required  a  highly  robust  density-energy  query  capability.  Our  mixing  rule  assumes  that  our 
materials  are  immiscible  and  are  in  mechanical  (pressure)  and  thermal  equilibrium.  In  this  mixing 
rule,  we  assume  that  we  have  a  mixture  of  species  equations  of  state  which  define  pressure  as  a 
function  of  density  and  temperature.  Using  the  immiscible  assumption,  we  find  the  volume  of  the 
mixture  is  the  sum  of  the  volume  occupied  by  each  species,  or 

i  NS  Y. 

-=  I  (3.5) 

P  f=l  P' 
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where  Y  is  the  species  mass  fraction  and  p,  is  the  density  of  pure  material  i  at  a  given 
temperature  and  pressure. 


p=pj(pi,T)Vie[\..NS], 


(3.6) 


where  p.(p,  ,7)  is  the  equation  of  state  for  species  i,  and  p  is  the  system  pressure.  The  energy  of 
mixture  is  given  by  the  mass  averaged  species  energies,  as  in 

NS 

e=  Z  Y.e.(T, pf),  (3.7) 

i=\ 


where  the  species  energy  equation  is  given  by 

ei (T , p] )  =  efj  + £ c,,, (t)c1t  +  [e, (t , p*  )-e;(7% J 


(3.8) 


where  the  last  term  of  this  expression  is  the  departure  function  that  accounts  for  the  effects  of 
non-ideality  in  the  equation  of  state.  Note,  equation  (3.8)  is  simply  a  trivial  regrouping  of  the 
energy  components  that  allows  us  to  view  the  energy  of  a  species  as  a  division  of  thermally 
perfect  and  thermally  imperfect  components.  Thus  the  departure  function  is  viewed  as  the 
deviation  of  the  energy  equation  from  a  thermally  perfect  gas.  The  departure  function  can  be 
computed  using  the  relation: 


[e/(r,A*)-^(7’L]=  1 


p-t 


dp(p'jj) 


dT 


dp, : 


Pi 


(3.9) 


Equations  (3.5),  (3.6),  and  (3.7)  describe  NS+2  nonlinear  equations  and  NS+  4  unknowns  (the 
unknowns  are  the  species  specific  volumes,  1/p*,  and  the  thermodynamic  variables  P ,  T ,  p ,  and 

e ).  Given  any  two  thermodynamic  state  variables,  this  system  of  equations  can  be  solved  to  find 
the  remaining  thermodynamic  variables  along  with  the  specific  volumes  implied  by  the  species 
pure  substance  densities,  p* .  Due  to  the  non-linearities  in  these  equations,  robustly  solving  for 

the  thermodynamic  state,  given  the  fluid  density  and  energy  (as  will  be  required  by  the  explicit 
time  integrator),  can  be  a  significant  challenge.  The  most  straightforward  approach  is  to  use  a 
multi-dimensional  Newton  method.  However,  when  the  initial  guess  is  far  away  from  the  final 
solution,  it  can  be  difficult  to  reliably  converge  to  a  solution.  We  have  developed  an  alternative 
approach  that  is  robust,  but  perhaps  with  a  sacrifice  in  computational  efficiency. 

To  obtain  a  robust  solution  to  the  above  non-linear  equations,  we  have  made  the  observation 
that  if  we  solve  for  pressure  and  temperature  instead  of  density  and  energy,  the  pressure  equation, 
i.e.,  equation  (3.6),  decouples  and  we  can  then  solve  NS  non-linear  scalar  equations.  Since  non¬ 
linear  root  finding  for  scalar  equations  can  be  made  robust  using  root  bracketing  techniques,  a 
pressure-temperature  query  has  a  very  robust  solution.  Can  we  utilize  this  observation  to  make  a 
robust  density-energy  EoS  query?  Yes,  through  an  indirect  route.  First,  we  can  obtain  a  robust 
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density-temperature  query  by  performing  a  bracketed-scalar  solve  for  the  pressure.  This  iterative 
solve  can  utilize  the  robust  pressure-temperature  query;  in  addition,  it  is  expected  that  pressure 
will  increase  monotonically  with  temperature.  Finally,  given  a  robust  density-temperature  query, 
a  robust  density-energy  query  can  be  formulated  by  solving  for  the  temperature  that  gives  the 
specified  energy.  All  of  these  queries  are  solved  using  a  bracketed  Newton  method  whereby  the 
values  that  bracket  a  solution  are  identified.  If  the  Newton  method  overshoots  the  bracket,  the 
robust  bisection  method  is  utilized  for  that  step.  The  resulting  EoS  query  evaluator  has  been 
shown  to  be  robust  in  practice. 

3.3  Species  equations  of  state 

This  section  will  outline  the  species  equations  of  state  implemented  for  blast  wave  modeling 
and  provide  a  description  of  the  definition  of  pressure  and  energy  for  each  species.  The  energy 
for  each  of  these  species  is  derived  using  the  departure  function  in  equation  (3.9). 

Perfect  gas 

The  perfect  gas  equation  of  state  provides  the  basic  equation  of  state  for  thermally  and  calorically 
perfect  materials.  The  equation  of  state  for  a  perfect  gas  is  given  by  the  relation: 

P(?,,T)=P,  ~T,  (3.10) 

i 

where  m.  is  the  molecular  mass  of  species  i  and  R=83\4.3J/(Kkmol )  is  the  universal  gas  constant. 
The  energy  is  given  by  the  relation: 


*  R 

e(pi,T)=n.  T+h, 

■>  J  ,  m  fi ’ 

i 


(3.11) 


where  hf]  is  the  heat  of  formation  for  species  i  and  is  defined  by  the  reference  enthalpy  at  the 

specified  reference  temperature.  For  air,  the  perfect  gas  model  is  specified  in  the  Loci/CHEM 
.  mdl  file  as 

//  Air  as  an  ideal  gas 

_Air  =  <  m=28.89,  n=2.5,  href=0,  sref=0,  Tref=298.0, 

Pref=10000.0  >  ; 


JWL  gas 

The  JWL  equation  of  state  defined  in  density-temperature  form  as  utilized  by  Lee  and  Tarver 
[Lee  and  Tarver,  1980]  defines  pressure  as  the  relation: 


/y.(p/,r)=dexp 


Pol 

f  P3 

r  0 

R\ -J 

+Bex p 

-R2-* 
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where  pQ  is  the  density  of  the  explosive  material  pre-detonation,  co,  A,  B,  R 1,  and  R2  are 
coefficients  of  the  JWL  EoS,  and  C  defines  the  gas  specific  heat  of  the  expanded  gas.  Note,  the 
default  units  for  A  and  B  are  Pascals,  while  C  is  Pascal/kelvin.  Utilizing  the  departure  equation, 
equation  (3.3.9),  the  energy  for  this  species  is  found  to  be  defined  by  the  following  relation: 
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A  species  is  prescribed  by  the  JWL  EoS  by  assigning  the  EoS  parameter  in  the  Loci/CHEM 
model  file.  For  example,  the  JWL  EoS  for  the  explosive  C-4  is  given  by  the  following  species 
model  definition: 

//  Explosive  gases  released  under  detonation  of  C-4 

_jwlgas  =  <m=45,n=3.3,  href=0 , sref=0 . 0 , Tref=0 , Pref=0 , mf=0  >  ; 

_jwlgas:  <eos=JWL (A=60 9 . 77e9 ,  B=12.95e9,  omega=0.25, 

Rl=4 . 5,R2=1 .4, 

rho_0=1601, Cv=le6) >  ; 

Note,  the  molecular  weight,  m  and  partition  function  parameter,  n,  do  not  need  to  be  set  to  a 
specific  values  as  they  will  be  redefined  to  match  the  co  and  Cv  provided  in  by  the  JWL  EoS. 

Elastic  solid 

The  equation  of  state  for  a  perfectly  elastic  solid  can  be  defined  by  the  simple  linear  barytropic 
equation  of  state  given  by  the  relation: 


P,{pi,T)=K 
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where  K  is  the  bulk  modulus  of  the  material,  and  p0  is  the  material  density  at  a  given  reference 
pressure,  P  .  Utilizing  equation  (3.3.9),  the  energy1  for  the  elastic  solid  EoS  is  given  by 

*  K-PO  K  *  R 

efS>i,T)-  *  +  ln(P  i)+nimT+hfi  (3-15) 

Pi  P0  i  J 

Note,  in  describing  an  elastic  solid,  the  following  properties  must  be  defined:  the  material 
density  at  a  given  reference  pressure,  the  material  bulk  modulus,  and  set  m  and  n  to  achieve  the 
material  C  .  For  example,  the  material  quartz  has  a  specific  heat,  C=133J/kgKm,  thus  if  we  set  the 


'Note:  We  cannot  integrate  to  a  zero  density  for  an  elastic  solid.  However,  the  relation  will  hold  for  any  arbitrary 
constant  of  integration  with  the  interpretation  that  we  are  integrating  from  some  given  reference  density,  p().  Since 

the  absolute  value  of  the  energy  is  not  critical  we  absorb  this  constant  of  into  the  heat  of  formation  term,  h  of  this 

Jl 

species. 
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molecular  weight  to  1  we  can  compute  n=CJR= 0.0882.  Since  quartz  has  a  density  of  2650/cg/m3 
and  a  bulk  modulus  of  £=36.4x1 09Pa,  the  model  for  quartz  using  an  elastic  solid  is  specified  as: 

/ /  Quartz  as  an  elastic  solid 

_quartz  =  <  m=l , n= . 0882 , href=0 ,  sref=0 . 0 , Tref=0 , Pref=0 , mf=0  >  ; 
_quartz  :  <eos=ElasticSolid (rho_0=2 650 , K=36 . 4e9 , P0=le5 )  >  ; 


Tait  EoS  for  liquid  water 

The  Tait  EoS  is  a  simple  and  accurate  barytropic  equation  of  state  for  liquid  water  under  high 
pressure.  This  equation  of  state  is  given  by  the  relation: 
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where  p()  is  the  material  density  at  reference  pressure  given  by  coefficient  A,  B  is  an  empirically 

determined  coefficient,  and  y  is  the  material  exponent.  This  equation  of  state  can  model  liquid 
water  with  the  following  coefficients:  A=\05Pa,  B=3.3\x\0HPa,  p(=l OOOkg/m3,  and  y=7.15.  The 

energy  as  defined  by  the  departure  function  for  this  EoS  is  expressed  by  the  relation: 
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For  this  equation  of  state,  the  material  specific  heat,  Cv  can  be  set  using  the  same  procedure  as 

described  for  the  elastic  solid.  The  model  input  to  describe  liquid  water  is  given  by  the  following 
material  definition: 

_water=  <m=l , n= . 50341 , href=0 , sref=0 , Tref=0 , Pref=0 , mf=0  >  ; 

_water:  <eos=Tait (A=le5 , B=3 . 31e8 , gamma=7 . 15 , rho_0=1000 )  >  ; 


3.4  Modeling  a  prescribed  explosive  burn 

We  have  developed  a  prescribed  bum  capability  that  can  be  used  to  simulate  the  propagation 
of  a  detonation  front  through  the  explosive  material.  The  current  model  assumes  that  the 
explosive  is  initiated  from  a  single  point  and  that  there  are  no  obstructions  to  the  detonation  front. 
In  the  prescribed  explosive  bum  methodology,  the  initiation  point  and  a  detonation  velocity  are 
provided  by  the  user.  The  explosive  bum  is  accomplished  by  converting  the  solid  explosive 
material  into  the  corresponding  gas  material  as  the  detonation  wave  passes  each  given  point  in 
the  mesh.  Usually  the  solid  explosive  before  detonation  is  modeled  as  an  elastic  solid,  while 
gasses  released  by  the  bum  are  modeled  using  the  JWL  EoS.  To  accomplish  the  appropriate 
energy  release  during  the  bum,  the  heat  of  formation  of  the  solid  explosive  material  is  set  such 
that  the  proper  heat  release  is  achieved. 
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The  detonation  front  is  computed  by  enforcing  a  burn  fraction  that  is  a  function  of  the  lighting 
time,  t  .  The  lighting  time  is  the  computed  time  that  the  detonation  front  will  arrive  at  a  given  cell 

and  is  computed  by  dividing  the  distance  to  the  initiation  point  by  the  detonation  velocity,  D.  The 
bum  fraction,  F,  is  zero  if  t<t{.  When  t>t ,  the  bum  fraction  is  defined  by 


F=max 


(3.18) 


where  t  is  the  current  simulation  time,  t  is  the  current  cell  lighting  time,  D  is  the  detonation 


velocity,  V=  is  the  relative  volume,  V 

P  ^ 


is  the  relative  volume  at  the  Chapman-Jouguet 


conditions,  and  A  is  an  estimated  grid  spacing  for  the  mesh  where  the  detonation  front  is 
propagating. 

The  explosive  bum  is  enabled  by  adding  the  explosiveBurnModel  option  to  the 
Loci/CHEM  .vars  file.  This  input  specifies  the  initiation  point,  the  detonation  velocity,  the 
Chapman-Jouguet  pressure,  and  the  species  that  describe  the  solid  and  gas  components  of  the 
explosive  material.  An  example  used  for  C-4  explosive  is  as  follows: 


/ /  initiate  in  center  of  charge 
explosiveBurnModel:  <  initiationPoint  = 
[-40mm, 0.0,0.005], D=8193m/s, P_CJ=28e9, 

solid=_c4 , gas=_j wlgas  > 


Creating  a  detonation  model 

In  order  to  determine  heat  release  from  the  solid  to  gas  species  decomposition,  the  correct 
heat  of  formation  must  be  determined  for  each  such  that  the  resulting  energy  release  is  achieved 
when  the  detonation  front  passes  a  particular  region..  Here  we  give  an  example  of  how  to  set  up 
these  parameters  for  the  material  C-4.  The  properties  form  the  JWL  EoS  for  C-4  are  as  follows: 
p0=160 lkg/m3,  D=S\93m/s,  PCJ=2Sx\09Pa,  E=9x\09J/m\  A=609J7x\09Pa,  B=l2.9x\09Pa, 

i?  1=4.5,  R2=\.4,  and  oo=0.25.  Note,  the  detonation  velocity,  D,  and  the  Chapman-Jouguet 
pressure,  P CJ>  are  provided  to  the  explosive  bum  input  given  in  the  .vars  file  as  specified 

above.  The  gas  species  is  covered  by  a  JWL  equation  of  state  (note  Cv  must  also  be  provided  to 

obtain  the  correct  specific  heat  of  the  expanded  detonation  products.)  The  example  is  repeated 
here: 


//  Explosive  gases  released  under  detonation  of  C-4 

_jwlgas  =  <m=45,n=3.3,  href=0 , sref=0 . 0 , Tref=0 , Pref=0 , mf=0  >  ; 

_jwlgas:  <eos=JWL (A=60 9 . 77e9 ,  B=12.95e9,  omega=0.25, 

Rl=4 . 5,R2=1 . 4, 

rho_0=1601, Cv=le6) >  ; 

Note  that  in  the  above  model,  setting  href=0  and  Tref=0  ensures  that  the  heat  of  formation 
for  this  species  is  identically  zero.  To  describe  the  explosive  bum,  we  must  also  specify  the  solid 
explosive  material  EoS.  For  this  we  utilize  the  elastic  solid  with  p0=1601,  A=10xl09,.P0=105, 

m=  1,  and  «=.18.  These  parameters  are  set  to  approximate  the  properties  of  the  solid  explosive 
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material  as  an  elastic  solid.  If  we  set  Tref =0  and  m=l,  then  we  can  compute  the  value  href  to 
release  E  units  of  energy  as  given  by  the  expression  (note  h^.  is  in  the  units  of  J/kmol  and  m  has 

the  units  kg/kmol): 


E-  (Kill  (p.) )£-/>, 


h=m 


-nRT 


Assuming  ro=300  Kelvin  we  can  then  compute  the  reference  enthalpy  at  absolute  zero 
because  the  solid  material  can  be  computed  as: 


9xl09-(l+ln(1601))10xl09-105 

h=  1 - 1 - V16Q1^ - -0.1 8(83 14.3)(300)=-47 159739  (3.20) 


Thus,  the  solid  C-4  material  will  be  defined  by  the  following  specification: 

c4=  <m=l,n=.18,  href=-47159739,  sref=0 . 0 , Tref=0 , Pref=0 , mf=0  >  ; 
c4 :  <eos=ElasticSolid (rho  0=1 601 , K=10e9 , P0=le5 )  >  ; 


3.5  Techniques  to  preserve  positive  mass  fractions 

When  simulating  multicomponent  flows,  it  is  possible  to  take  a  time-step  such  that  the 
material  in  a  cell  is  completely  depleted,  resulting  in  the  time  evolution  of  negative  mass 
fractions  of  material.  This  non-physical  circumstance  is  unacceptable  because  advancing  the 
solution  becomes  nearly  impossible.  For  first-order  spatial  approximations,  depletion  of  the 
material  from  a  cell  is  avoided  if  a  time-step  is  employed  that  satisfies  the  stable  CFL  condition. 
However,  when  reconstructing  higher-order  mixture  fractions,  it  is  possible  to  deplete  a  material 
even  when  using  a  much  smaller  time-step  To  allow  for  second-order  spatial  reconstruction  of 
species  mass  fractions,  we  first  limit  the  mass  fraction  extrapolations  so  that  a  negative  mass 
fraction  will  not  occur  (assuming  a  first-order  upwind  convection).  While  such  limiting  helps 
avoid  the  evolution  of  negative  mass  fractions,  they  can  still  evolve.  Therefore,  to  provide  an 
effective  strategy  that  avoids  negative  mass  fraction  evolution  in  the  time-stepping  algorithm,  we 
limit  the  stable  time-step  to  include  the  time  required  to  deplete  a  cell  of  all  species  material 
using  the  first  time-step  residual..  However,  for  multicomponent  simulations,  a  CFL  setting 
between  0.75  and  0.5  still  may  be  required  to  avoid  the  evolution  of  negative  mass  fractions. 

3.6  Validation  studies 

Hemispherical  surface  explosion 

To  test  our  JWL  EoS  and  our  explosive  bum  model,  we  compared  the  results  of  simulations 
of  a  rigid,  surface-laid  hemispherical  charge  of  TNT  with  a  radius  of  144mm  and  a  mass  of 
10.19kg,  and  we  compared  the  blast  strength  and  impulse  with  results  provided  in  the  Fiserova 
PhD  dissertation  [Fiserova,  2006],  The  resulting  blast  wave  parameters,  namely  overpressure  and 
specific  impulse,  were  examined  at  distances  between  200 mm  and  700 mm  above  the  ground. 
JWL  equation  of  state  and  perfect  gas  equation  of  state  for  modeling  the  explosive  gas  were  both 
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evaluated.  Two  different  cases  using  the  JWL  EoS  were  studied:  one  utilizing  the  prescribed 
explosive  bum  detonation  model  and  the  other  utilizing  a  bursting  high  pressure  sphere  with 
initial  pressure  given  by  the  equilibrium  pressure  of  the  volume-confined  explosive  gases.  Our 
simulation  findings  were  compared  with  empirical  data  obtained  from  CONWEP  and  simulation 
results  from  AUTODYN  provided  by  Fiserova  [Fiserova,  2006]. 

When  the  explosive  bum  model  was  used  in  this  study,  the  TNT  solid  material  was  modeled 
as  an  elastic  solid  with  a  bulk  modulus  K=\0GPa  and  density  at  ambient  conditions  of 
p0=1630 kg/m3.  The  explosive  gas  was  modeled  using  the  JWL  equation  of  state  with 

A=31\.2GPa,  B=3.23lGPa,  oo=0.3,  Rl=4.15,  and  R2=0.95.  Air  was  modeled  as  a  perfect  gas 
with  a  y=1.4  and  an  average  molecular  weight  of  m= 28.89.  The  explosive  bum  was  initiated  from 
the  center  of  the  sphere  that  formed  the  hemispherical  charge  with  a  detonation  velocity  of 
D=6900m/s  and  a  Chapman- Jouguet  pressure  of  P  =2 1  GPa. 

Grid  convergence  was  studied  using  5 mm,  3mm,  1  mm,  0.15mm  and  0.315mm  cell  sizes. 
Figures  3.1,  3.2  and  3.3  are  time  histories  of  the  overpressure  obtained  from  the  JWL  EoS  with 
detonation  model,  the  JWL  EoS  without  detonation  model,  and  the  perfect  gas  EoS,  respectively. 
The  overpressure-time  histories  were  recorded  at  distances  between  200mm  and  100mm  above 
the  ground  at  100mm  increments.  For  the  models  using  the  perfect  gas  EoS  and  the  JWL  EoS 
without  detonation,  it  was  observed  that  the  cell  size  had  little  effect  on  the  resulting  blast  wave 
parameters  except  at  the  closest  distance  above  the  ground.  However,  the  results  using  the  JWL 
EoS  with  the  explosive  bum  model  were  shown  to  be  more  dependent  upon  the  grid  size.  For  cell 
sizes  of  0.5mm  and  0.375mm,  the  maximum  over-pressures  are  convergent.  For  the  finest  grid 
size  in  JWL  EoS  with  detonation,  oscillations  in  the  pressure  after  the  initial  shock  are  apparent 
in  the  results.  Investigation  of  these  oscillations  showed  that  they  were  caused  by  apparent 
Richtmyer-Meshkov  and  Rayleigh-Taylor  instabilities  forming  in  the  contact  discontinuity 
between  the  explosive  products  and  air  as  is  shown  in  Figure  3.4.  Some  of  these  instabilities  seem 
to  have  been  seeded  by  irregularities  created  by  the  unstmctured  mesh  used  in  the  simulation. 
However,  such  turbulent  mixing  of  the  contact  discontinuity  is  physical  and  expected  based  on 
high  speed  photographic  recordings  of  the  fireball  in  explosive  events. 


Figure  3.1:  Time  history  of  overpressure  for  TNT  charge  with  explosive  burn  model  using 

JWL  EOS. 
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Figure  3.2:  Time  history  of  overpressure  for  TNT  charge  without  explosive  burn  model 

using  JWL  EOS. 


Figure  3.3:  Time  history  of  overpressure  using  perfect  gas  EOS. 
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Figure  3.4:  Simulated  instabilities  fireball  on  the  0.375mm  mesh. 

Blast  wave  parameters  obtained  from  our  simulation  results  were  compared  with  empirical 
data  from  CONWEP  and  simulation  findings  from  AUTODYN,  as  shown  in  Figs  3.5  and  3.6. 
For  the  perfect  gas  EoS  and  the  JWL  EoS  without  detonation  model,  results  from  the  1  mm  grid 
size  were  utilized.  For  the  JWL  EoS  with  detonation  model,  results  from  the  0.5 mm  grid  size 
were  used.  For  the  purpose  of  fair  comparisons,  results  from  the  0.5 mm  grid  size  were  chosen  for 
the  AUTODYN  simulations.  Figure  3.5  presents  the  maximum  over-pressures  at  different 
distances  above  the  ground.  At  the  closest  point  above  the  ground  (200 mm),  all  the  simulation 
results  over-predict  by  at  least  37.5  percent.  However,  it  is  noted  that  CONWEP  is  unable  to 
provide  data  for  distances  below  300 mm  in  this  case;  therefore,  the  blast  parameters  from 
CONWEP  at  distances  below  300 mm  are  evaluated  by  extrapolation.  At  distances  above  300 mm, 
the  simulations  provide  reasonable  results  with  the  perfect  gas  EoS  and  the  JWL  EoS  without 
detonation  underestimating  the  maximum  over-pressure,  and  the  JWL  EoS  with  detonation  and 
the  AUTODYN  overestimating  it.  Figure  3.6  shows  the  specific  impulse  at  different  distances 
above  the  ground  using  different  models.  The  specific  impulse  is  defined  as  the  time  integration 
of  the  overpressure.  Again,  the  simulation  results  from  the  models  over-predict  when  compared 
with  the  CONWEP  extrapolated  data  at  the  closest  distance.  At  the  farthest  distances  (larger  than 
500mm),  the  simulation  findings  appear  to  be  lower  than  experimental  data. 


44 


UNCLASSIFIED 


Figure  3.5:  Comparison  of  maximum  overpressure  using  different  models. 


Figure  3.6:  Comparison  of  specific  impulse  using  different  models. 


Small-scale  explosion  in  dry  sand 

The  second  validation  study  simulated  the  explosion  of  C-4  charges  buried  under  dry  sand. 
Three  burial  depths,  were  given  in  terms  of  the  thickness  of  soil  overlaying  the  mine:  0 mm  for 
flush  burial,  and  3 cm  or  8 cm  overburdens  for  buried  landmines.  Blast  wave  parameters  were 
measured  in  the  computations  at  stations  located  30cm  and  70 cm  above  the  ground.  The 
explosive  charge  was  defined  by  100  grams  of  C-4  in  a  disk-shaped  charge  with  a  radius  of  3.2 cm 
and  a  thickness  of  2 cm.  The  explosive  charge  was  initiated  in  the  center.  The  solid,  unexploded 
C-4  was  modeled  as  an  elastic  solid  with  a  bulk  modulus  of  K=\0GPa  and  a  density  at  ambient 
conditions  of  p0=1601  kg/mi.  The  explosive  gasses  were  modeled  using  the  JWL  EoS  with 

A=609 .77 GPa,  B=\2.95GPa,  oo=0.25,  7?1=4.5,  and  R2=\.4.  The  explosive  bum  was  modeled 
with  a  detonation  velocity  D=8 1 93 m/s  and  a  Chapman-Jouguet  pressure  P (,  =2%GP a.  Dry  sand 

was  modeled  as  a  mixture  of  solid  quartz  and  air,  with  solid  quartz  being  modeled  as  an  elastic 
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solid  with  a  bulk  modulus  of  K=33.11GPa  and  a  density  at  ambient  conditions  of  p0=2650 kglmi. 

The  mixture  contains  .0311%  air  by  mass.  The  resulting  mixture  density  was  the  same  at  ambient 
conditions  as  the  sand  utilized  in  the  experiment.  Notably,  in  this  simulation  there  were  four 
species:  air,  quartz,  jwlgas,  and  C4.  Sand  was  created  by  mixing  quartz  and  air  in  appropriate 
proportions.  As  the  detonation  front  expanded,  the  C4  species  converted  to  jwlgas. 

Simulations  of  the  buried  charge  were  carried  out  on  three  grids,  starting  with  one  that  had  a 
mesh  spacing  of  1  mm  in  the  explosive  charge,  and  then  successively  refining  said  grid  to  achieve 
0.5 mm  and  0.25 mm  mesh  spacings.  Pressure  probes  were  placed  at  30cm  and  70 cm  above  the 
ground.  At  the  time  of  reporting,  only  the  flush  cases  had  been  run  on  all  three  meshes  for  a 
sufficient  amount  of  time  to  determine  the  impulse  at  the  70mm  probe.  For  the  other  depths  of 
burial  (DOB),  the  simulations  should  be  run  longer  to  get  a  complete  grid  convergence 
comparison.  Results  from  the  flush  case  and  results  from  the  30cm  probes  on  the  other  cases 
indicate  that  the  1mm  mesh  is  satisfactorily  converged.  Comparison  with  experiments  and  other 
codes  will  therefore  be  given  for  the  results  of  simulations  on  the  1mm  mesh. 

Results  from  the  1mm  mesh  for  the  mine  that  was  buried  flush  with  the  ground  for  times  50p.v, 
lOOps,  500p^,  lOOOps,  1500p^,  and  2000|xs'  are  given  in  Figures  3.7-3.12,  respectively.  These 
figures  show  an  upward  directed  blast  wave  with  significant  compaction  of  the  sand  underburden 
in  the  early  stages  of  the  blast  formation.  As  the  explosive  blast  evolves,  the  crater  forms  an 
upward  jet  of  exploding  gases.  These  gases,  which  expand  behind  the  leading  shock  wave, 
become  unstable  and  mix  in  a  turbulent  fashion  after  t=500p.v. 

Results  from  the  simulation  on  the  1mm  mesh  for  the  mine  that  was  buried  under  3 cm  of  sand 
for  time  50p5,  lOOps,  500p5,  lOOOps,  1500^5,  and  2000ps  are  given  in  Figures  3.13-3.18, 
respectively.  These  figures  show  the  compaction  and  deformation  of  the  surface  sand  in  response 
to  the  growing  bubble  of  explosive  gases.  By  t=1000ps  the  sand  layer  begins  to  become  unstable 
and  mix  with  the  explosive  gases,  and  by  t=2000|xs  the  sand  interface  breaks  into  disconnected 
regions,  as  can  be  observed  in  Figure  3.18. 

Results  from  the  simulation  on  the  1mm  mesh  for  the  mine  that  was  buried  under  8 cm  of  sand 
for  time  50ps,  lOOps,  500ps,  lOOOps,  1500^,  2000ps,  3000p^,  and  4000p^  are  given  in  Figures 
3.19-3.26,  respectively.  At  this  DOB,  the  simulations  show  a  significant  bubble  of  explosive 
gases  forming  under  the  sand  before  eruption  through  the  surface.  By  /=2000pv,  instabilities  in 
the  outer  sand  layer  are  beginning  to  cause  mixing  between  the  sand  and  the  expanding  explosive 
gases.  By  /=4000pv,  the  explosive  gases  have  broken  through  the  sand  layer  and  there  is 
significant  mixing  evident  in  the  fireball.  All  of  these  simulations  are  qualitatively  similar  to 
simulations  documented  by  Fiserova  [Fiserova,  2006],  Wang  [Wang  2001],  and  Grujicic 
[Grujicic  et  al.,  2006]. 

We  can  show  results  on  the  0.25mm  spacing  grid  that  illustrate  some  of  the  fine  details 
captured  in  our  simulations.  Figure  3.27  shows  the  density  plot  from  this  high  resolution 
simulation  at  1400p,v  after  detonation.  Here  we  see  that  the  overburden  sand  that  is  thrust  in  the 
air  by  the  blast  undergoes  severe  mixing,  which  suggests  the  presence  of  large  instabilities.  These 
instabilities  appear  to  come  from  baroclinic  torque  caused  by  pressure  disturbances  moving 
through  the  thin,  high-density  region  that  forms  from  the  soil  overburden.  In  turn,  these  rapidly 
rotating  regions  produce  some  rather  violent  motion  of  the  sand  ejecta  that  results  in  regions  of 
high  and  low  pressure  and  which  further  enhance  the  unstable  motion.  These  effects  can  be  seen 
even  more  dramatically  in  the  numerical  Schlieren  shown  in  Figure  3.28.  We  suspect  that  this  is 
not  physically  real,  as  the  sand  should  become  a  dispersed  phase  with  momentum  somewhat 
decoupled  from  the  gas  phase  pressure;  however,  the  current  homogeneous  multiphase  mixture 
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does  not  permit  this  to  occur.  It  is  suggested  that  removing  sand  that  has  become  sufficiently 
dispersed  and  placing  it  in  a  Lagrangian  multiphase  model  would  improve  the  physical  realism  of 
the  model  at  this  stage  of  the  computations.  However,  given  sufficient  mesh  resolution,  it  appears 
that  the  sand  would  naturally  break  into  a  dispersed  material  cloud  due  to  the  unstable 
mechanisms  that  we  observe  in  these  simulations. 


t=50us 


Figure  3.7:  Density  plot,  DOB=Oc/w,  t=50ps. 
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Figure  3.8:  Density  plot,  DOB=0 c/w,  t=100|xs. 
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Figure  3.9:  Density  plot,  DOB=Ocam,  t=500ps. 


Figure  3.10:  Density  plot,  DOB=Ocam,  t=1000ps. 
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Figure  3.11:  Density  plot,  DOB=Ocam,  t=1500ps. 


Figure  3.12:  Density  plot,  DOB=Ocam,  t=2000ps. 
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Figure  3.13:  Density  plot,  DOB=3 cm,  t=50(j.s. 
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Figure  3.14:  Density  plot,  DOB=3 c/w,  t=100|is. 
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Figure  3.15:  Density  plot,  DOB=3c/w,  t=500|xs. 


Figure  3.16:  Density  plot,  DOB=3cam,  t=1000ps. 
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Figure  3.17:  Density  plot,  DOB=3c/«,  t=1500ps. 


Figure  3.18:  Density  plot,  DOB=3cm,  t=2000(j.s. 
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Figure  3.19:  Density  plot,  DOB=8cm,  t=50(j.s. 


Figure  3.20:  Density  plot,  DOB=8cm,  t=100ps. 
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t=500us 


density  (tg/m'M) 


Figure  3.21:  Density  plot,  DOB=8cm,  t=500|xs. 


Figure  3.22:  Density  plot,  DOB=8c/w,  t=1000|Ls. 
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Figure  3.23:  Density  plot,  DOB=8cm,  t=1500ps. 


Figure  3.24:  Density  plot,  DOB=8c/w,  t=2000ps. 
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Figure  3.25:  Density  plot,  DOB=8cm,  t=3000ps. 


Figure  3.26:  Density  plot,  DOB=8c/w,  t=4000ps. 
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Figure  3.27:  Density  plot  for  dry  sand,  DOB=3cn/,  t=1400ps,  grid  spacing=0.25«zAM. 


Figure  3.28:  Numerical  Schlieren,  DOB=3cam,  t=1400ps,  grid  spacing=0.25n/AM. 
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In  addition  to  these  qualitative  assessments,  a  case  with  saturated  sand  was  also  run.  For 
saturated  sand,  the  air  volume  was  replaced  with  water,  modeled  using  the  Tait  EoS;  however,  we 
included  a  very  small  fraction  of  air  to  allow  for  material  expansion  because  the  Tait  EoS  does 
not  properly  include  cavitation  effects.  Thus,  the  saturated  sand  model  consists  of  a  quartz  mass 
fraction  of  0.77,  a  water  mass  fraction  of  0.229995,  and  an  air  mass  fraction  of  5xl0-6.  For 
comparison,  the  same  buried  mine  simulation  for  a  DOB  of  3 cm  was  run  for  both  dry  sand  and 
saturated  sand.  The  pressure  plots  for  these  simulations  at  t=500p5  are  given  in  Figures  3.29  and 
30.  Note  in  these  figures,  the  white  curve  indicates  the  sand  interface  in  order  to  make  the 
location  of  the  sand  clear.  For  the  dry  sand  in  Figure  3.29,  the  shock  moves  slowly  due  to  the  low 
sound  speed  in  dry  sand.  The  detonation  products  are  less  focused  in  an  upward  direction  than  the 
saturated  sand  case  shown  in  Figure  3.30.  In  the  saturated  sand  case,  a  fast  moving  shock  wave  is 
evident  in  the  sand,  reflecting  the  fact  that  saturated  sand  has  a  much  higher  sound  speed.  In 
addition,  the  expanding  detonation  products  are  more  focused  in  an  upward  direction  than  the  dry 
sand  case.  Qualitatively,  these  results  are  very  similar  to  the  simulations  performed  by  Grujicic 
[Grujicic  et  al.,  2006], 

For  quantitative  validations,  the  simulation  results  from  the  1  mm  grid  are  compared  with 
experimental  data  obtained  by  Bergeron  et  al.  [Bergeron  et  al.,  1998],  In  addition,  the 
AUTODYN  results  from  Fiserova  [Fiserova,  2006]  and  LS-DYNA  results  from  Wang  [Wang 
2001]  are  shown  for  comparison.  While  a  final  grid  convergence  study  was  not  fully  complete  at 
the  time  this  report  was  generated,  preliminary  results  indicate  that  the  1  mm  meshes  were  of 
sufficient  resolution  to  accurately  model  the  problem  at  hand.  For  the  quantitative  comparison, 
the  following  blast  wave  parameters  were  considered:  time  of  arrival  of  the  blast  wave  front, 
maximum  overpressure,  and  specific  impulse.  The  results  from  experiments  and  the  simulations 
for  flush  3 cm  DOB  are  presented  in  tables  3.1  through  3.6.  For  flush  deployment,  the  simulations 
results  were  similar  to  the  AUTODYN  results  and  also  in  reasonable  agreement  with  the 
experimental  data  in  terms  of  impulse  and  arrival  time.  In  terms  of  overpressure  predictions,  the 
simulations  over-predicted  by  an  amount  similar  to  AUTODYN,  particularly  for  the  close  probe 
located  30cm  above  the  ground.  For  the  30cm  DOB,  the  simulations  again  produced  results 
similar  to  AUTODYN.  The  simulations  matched  the  experimental  results  reasonably  well  in 
terms  of  all  parameters:  maximum  overpressure,  specific  impulse,  and  arrival  time. 
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Figure  3.29:  Pressure  plot  for  dry  sand,  DOB=3cm,  t=500ps. 


Figure  3.30:  Pressure  plot  for  saturated  sand,  DOB=3c»i,  t=500ps. 


Table  3.1:  Comparison  of  maximum  overpressure  for  flush  deployment  in  dry  sand  for 

lOOg  C-4  charge. 


Position 

flush,  DOB=0  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

2797 

8160 

7380 

1359 

700 

1189 

1583 

1409 

580.8 
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Table  3.2:  Comparison  of  specific  impulse  for  flush  deployment  in  dry  sand  for  lOOg  C-4 

charge. 


Position 

Flush,  DOB=0  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

85.8 

98.0 

98.0 

86 

700 

116.4 

177 

169.6 

137.5 

Table  3.3:  Comparison  of  arrival  time  for  flush  deployment  in  dry  sand  for  lOOg  C-4 

charge. 


Position 

Flush,  DOB=0  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

94.8 

71.3 

76 

90 

700 

285.6 

271.6 

295 

300 

Table  3.4:  Comparison  of  maximum  overpressure  for  3 cm  buried  deployment  in  dry  sand 

for  lOOg  C-4  charge. 


Position 

Buried,  DOB=30  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

724.8 

845 

929.7 

613.3 

700 

304.5 

462.3 

334.1 

290.1 

Table  3.5:  Comparison  of  specific  impulse  for  3 cm  buried  deployment  in  dry  sand  for  lOOg 

C-4  charge. 


Position 

Buried,  DOB=30  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

106.8 

90.7 

86.1 

174.5 

700 

57.2 

43.1 

68.0 

77.9 

Table  3.6:  Comparison  of  arrival  time  for  3 cm  buried  deployment  in  dry  sand  for  lOOg  C-4 

charge. 


Position 

Buried,  DOB=30  mm 

above  ground 
(mm) 

Measured 

Our  results 

AUTODYN 

LS-DYNA 

300 

266 

356 

318 

270 

700 

784 

804 

740 

710 
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Figures  3.31-3.37  show  comparisons  of  displacement  parameters  between  our  simulation 
results  and  the  experimental  data  obtained  by  Bergeron  et  al.  [Bergeron  et  al.,  1998],  In  addition, 
the  AUTODYN  results  from  Fiserova  [Fiserova,  2006]  and  LS-DYNA  results  from  Wang  [Wang 
2001]  are  shown  for  comparison.  The  explosive  charge  is  lOOg  C-4  buried  under  dry  sand.  The 
availability  of  experimental  data  to  compare  dictated  running  two  DOB  cases,  namely 
DOB=30mm  and  DOB=Omm.  From  the  grid  convergence  study  run  with  1  mm,  0.5 mm  and 
0.25 mm  grid  sizes,  we  found  that  grid  convergence  is  reached  with  0.5 mm  grid  size.  Therefore, 
the  results  from  the  0.5 mm  grid  size  were  utilized  in  all  comparisons.  Vertical  displacement  of 
the  soil,  which  is  measured  as  the  distance  between  the  ground  and  the  ejecta  front,  is  plotted  in 
Figure  3.31.  Predictions  from  our  simulations  agree  well  with  the  experimental  data.  Figure  3.32 
shows  the  temporal  variation  in  the  crater  diameter.  Fairly  good  agreement  is  achieved  between 
the  simulation  and  the  measurement.  At  later  times,  our  simulations  over-predict  the  crater 
diameter  by  approximately  30  percent.  It  is  likely  that  at  later  times,  the  lack  of  a  soil  strength 
model  is  responsible  for  this  discrepancy;  specifically,  the  strain  rates  diminish,  which  increases 
the  relative  importance  of  the  material  strength  in  the  material  kinematics.  Figures  3.33  and  3.34 
show  the  height  and  width  of  the  detonation  products  cloud,  respectively.  The  agreement  between 
simulation  results  and  experimental  data  is  very  good  in  both  cases.  We  note  that  at  longer  times, 
our  model  under-predicts  values  relative  to  the  experimental  data.  We  suspect  that  this  may  be 
due  to  a  lack  of  a  failure  mechanism  for  the  soil.  In  our  simulations,  the  soil  encloses  the  growing 
cloud  of  detonation  products,  whereas  in  the  AUTODYN  simulations  the  soil  ruptures  thereby 
releasing  the  explosive  products.  In  our  method,  the  soil  remains  in  velocity  equilibrium  with  the 
explosive  gases.  Rather  than  failure,  we  propose  that  the  soil  should  transition  into  a  multiphase 
cloud  and  allow  the  solid  soil  material  to  take  on  independent  velocity  groups  from  the  gas  phase. 
We  are  now  extending  our  model  to  include  a  multiphase  modeling  capability  where,  instead  of 
failure  removing  material,  failure  would  inject  multiphase  particles  into  the  simulation.  Such  a 
model  would  allow  for  decoupling  between  the  soil  particles  and  gas  velocities,  thereby 
presumably  creating  a  more  rapid  growth  of  the  detonation  cloud. 


Figure  3.31:  Time  history  of  the  height  of  ejecta  at  DOB=30mm. 
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Figure  3.32:  Time  history  of  the  width  of  crater  at  DOB=30mm. 


Figure  3.33:  Time  history  of  the  height  of  the  detonation  product  at  DOB=30mm. 
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Figure  3.34:  Time  history  of  the  width  of  the  detonation  product  at  DOB=30mm. 

Figures  3.35-3.37  show  the  results  for  the  case  of  zero  DOB  (the  top  of  the  mine  is  at  ground 
level).  Figure  3.35  shows  the  crater  diameter,  which  agrees  with  experimental  results  very  well. 
Figure  3.36  and  3.37  show  the  height  and  width  of  the  detonation  products  cloud,  respectively. 
The  trends  of  our  simulation  agree  with  those  from  the  experiment.  We  note  that  in  the  case  of  a 
flush  buried  mine,  there  would  be  little  sensitivity  to  failure  because  there  is  no  soil  covering  the 
mine  to  rupture.  In  addition,  the  crater  formation  will  be  mostly  dominated  by  soil  compression 
rather  than  shear,  so  we  would  expect  little  sensitivity  in  these  results  to  material  strength.  In 
addition,  since  our  algorithms  employ  numerical  methods  that  conserve  total  energy  (including 
material  kinetic  energy),  we  will  capture  the  gas  dynamics  very  well,  even  on  a  reasonably  coarse 
mesh.  Thus,  in  many  ways  the  flush  case  is  optimal  for  a  simulation  using  numerical  methods 
that  are  optimized  for  gas  dynamics  such  as  those  employed  in  our  models. 


Figure  3.35:  Time  history  of  the  width  of  crater  at  DOB=Omm. 
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Figure  3.36:  Time  history  of  the  height  of  the  detonation  product  at  DOB=Omm. 


Figure  3.37:  Time  history  of  the  width  of  the  detonation  product  at  DOB=Omm. 


Mine  impulse  pendulum  (MIP)  explosion  in  prairie  soil 

The  mine  impulse  pendulum  is  a  horizontal  ballistic  pendulum  that  measures  the  effective 
impulse  from  a  landmine  explosion.  The  experiment  was  performed  using  DRDC  Suffield’s 
small  pendulum  apparatus  [Coffey  et  ah,  1998],  The  pendulum  consists  of  a  target  plate  attached 
to  two  rotating  steel  beams  affixed  to  a  stationary  base.  The  pendulum  target  plate  is  suspended  at 
a  specific  standoff  distance  over  a  buried  landmine.  The  landmine  is  detonated;  maximum  angle 
of  inclination  obtained  by  the  pendulum  arm  is  measured;  and  effective  mechanical  impulse 
applied  to  the  target  by  the  landmine  can  be  derived  as  a  function  of  the  maximum  inclination 
angle  of  the  pendulum.  Details  describing  the  MIP  experiment  can  be  found  in  the  PhD 
dissertation  of  Fiserova  [Fiserova,  2006], 
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In  our  simulation,  an  axisymmetric  plate  is  used  to  represent  the  square  witness  plate  in  the 
experimental  setup.  The  axisymmetric  disk  has  the  surface  area  and  moment  of  inertia  equivalent 
to  those  of  the  square  target  plate  in  the  experiment.  The  dimension  of  the  disk  is  137 6mm  in 
diameter  and  68.76 mm  in  height,  with  standoff  distance  of  400/h/h  above  the  ground.  An 
explosive  charge  of  1  Kg  C-4  is  buried  with  one  of  two  overburdens  (DOB=  \ 00mm,  50 mm).  The 
diameter  and  thickness  of  the  explosive  charge  are  134 mm  and  45 mm,  respectively.  The  soil 
properties  of  the  trials  used  for  validations  in  this  section  are  shown  in  Table  3.7.  The  density 
values  for  air,  water  and  solid  in  the  soil  are  \2kglnr’,  1000 kg/m3  and  2650 kg/m3,  respectively. 
The  value  of  dry  density  of  the  soil  is  set  to  1322  kglm 3  for  all  the  trials  in  our  simulations.  Using 
the  wet  and  dry  density  values,  the  mass  fractions  of  the  three  phases  (air,  water,  solid) 
composing  the  soil  can  be  determined,  which  are  shown  in  Table  3.8,  where  mfs,  mf ,  and  mfa 


represent  the  mass  fraction  of  solid,  water  and  air  in  the  soil,  respectively.  The  water  content 
percentage  was  determined  experimentally  by  measuring  the  mass  of  the  soil  and  then  drying  the 
soil  and  measuring  the  mass  of  the  dry  soil.  The  moisture  content  is  then  defined  as 


w  =  ■ 


(muet  ~™dry) 


m 


(3.21) 


dry 


The  mass  fractions  computed  in  Table  3.8  were  derived  by  assuming  that  the  void  volume  is 
conserved  and  that  the  dry  density  and  wet  density  can  be  written  as 


Pdry  ~ 


P solid  V solid  P  air  ^ air 

V  +v 

solid  air 


’  P  wet 


P solid  ^ solid  P  water  ^ water  P air  ^ air 

v  +  v  +v 

solid  water  air 


(3.22) 


The  mass  fractions  were  then  derived  by  considering  (while  neglecting  the  air  contributions)  that 
PWet  =  Pdry  (1  +  w)  •  Combining  these  equations  and  selecting  an  arbitrary  reference  for  Vsolid ,  the 


remaining  volume  fractions  were  computed  by  assuming  that  the  pore  volume  is  conserved: 
[Kir  L ry  =  \Fwater  +  Kir  \Wet  •  From  these  volumes,  component  mass  fractions  were  computed. 


For  the  present  simulation,  the  solids  were  given  the  material  properties  of  quartz.  Since  this 
type  of  soil  contains  sand  and  clay  particles,  this  assumption  is  questionable.  However,  based  on 
presently  available  data,  the  material  properties  of  isolated  clay  (without  the  voids  and  water)  are 
unknown.  The  explosive  gas  was  described  by  the  JWL  equation  of  state,  and  the  air  was 
assumed  to  be  an  ideal  gas.  Impulse  exerted  on  the  target  plate  was  computed  and  compared  with 
the  experimental  data. 


Table  3.7:  Soil  properties  for  simulated  trials. 


Moisture  Content  (%) 

Wet  density  (kg/m3) 

6 

1401 

7.7 

1423 

12.2 

1483 

15.0 

1520 

18.2 

1562 

20.9 

1598 

23.7 

1635 

28.7 

1704 

65 


UNCLASSIFIED 


Table  3.8:  Calculated  mass  fraction  for  simulated  trials. 


Moisture  content  (%) 

mfs 

mf 

mf 

6 

0.9432598 

0.05669 

0.0003789 

7.7 

0.928165 

0.0715867 

0.0003304 

12.2 

0.891089082 

0.108892255 

0.000268851349 

15.0 

0.8694 

0.130628 

0.00023305 

18.2 

0.845755041 

0.154181466 

0.000194178589 

20.9 

0.826996 

0.173127 

0.0001630328 

23.7 

0.80804956 

0.191823766 

0.0001321024 

28.7 

0.7768 

0.2233 

8.0307298e-05 

Figure  3.38  shows  a  comparison  of  impulse  on  the  target  plate  versus  soil  moisture  content  for  a 
50 mm  overburden  landmine  explosion.  Simulation  results  from  AUTODYN  and  Martec’s 
Chinook  code  [Martin  and  Link,  2003]  are  also  included  for  comparison  purposes.  We  note  that 
the  formulation  of  the  Chinook  code,  which  is  similar  to  our  methodology,  employs  algorithms 
typically  used  to  solve  inviscid  gas  dynamics  equations.  The  figure  demonstrates  that  the  trends 
of  increasing  impulse  as  a  function  of  moisture  content  are  reproduced  and  qualitatively 
consistent  with  the  experimental  data.  However,  our  simulations  over-predict  the  impulse  by  5  to 
50  percent  in  comparison  with  the  experimental  data.  Results  from  AUTODYN  are  smaller  in 
magnitude  than  our  predictions  and,  in  general,  show  good  agreement  with  the  experimental  data. 
Similar  observations  can  be  made  for  the  100 mm  overburden  case,  as  depicted  in  Figure  3.39. 

The  over-prediction  of  the  impulse  for  our  present  model  is  not  unexpected,  as  we  will  explain 
shortly.  To  understand  the  details  of  this  prediction,  we  studied  the  evolution  of  the  plate  loading 
for  the  18.2%  moisture  content  and  50 mm  overburden  case.  Figure  3.40  shows  the  integrated 
momentum  flux  on  the  plate  as  a  function  of  time.  We  note  that  there  are  three  spikes  in 
momentum  transfer:  the  first  at  t=0.5ms,  the  second  at  t=0.9ms,  and  the  third  at  t=\.lms.  Between 
t=0.5ms  and  t=0.9ms,  soil  impacts  over  an  increasingly  large  area,  combined  with  a  reduced 
momentum  transfer  in  the  center  of  the  plate  due  to  a  soil  rebound  effect.  By  t=0.9ms,  as 
illustrated  in  Figure  3.41,  the  explosive  gases  are  enclosed  by  the  soil.  The  soil  velocity  at  this 
time  is  primarily  in  the  direction  normal  to  the  plate.  The  dip  in  the  momentum  transfer  observed 
at  t=\.2ms  is  due  to  the  rebound  of  the  soil  that  struck  the  plate  at  t=0.9ms,  while  the  final  peak  at 
t=\.lms  is  caused  by  an  increasingly  glancing  blow  of  the  soil  and  gas  material  with  the  target 
plate,  as  can  be  seen  in  Figure  3.42.  Also  noteworthy  is  the  presence  of  a  thin  region  of  dense, 
high-pressure  soil  material  in  the  region  immediately  adjacent  to  the  witness  plate  throughout  the 
duration  of  peak  momentum  transfer  to  the  plate.  Although  in  the  present  model  we  cannot 
distinguish  between  the  contributions  to  the  loading  due  to  the  soil  and  the  gases,  the  simulations 
support  that  a  significant  proportion  of  the  loading  is,  in  fact,  due  to  the  large  density  (and 
consequently  large  momentum)  of  the  soil  material.  That  our  model  is  a  homogeneous  mixture  of 
quartz,  air,  and  water  implies  that  the  soil  must  have  the  same  velocity  as  the  gas  phase.  As  a 
result,  the  model  transfers  too  much  momentum  to  the  soil.  A  better  model  for  this  case  would 
allow  the  soil  to  become  a  multiphase  cloud  of  particles  at  an  appropriate  time.  Once  the  soil 
particles  are  free  to  have  their  own  velocity,  the  momentum  coupling  will  become  less  efficient. 
Therefore,  there  will  be  less  momentum  in  the  soil  that  strikes  the  plate,  reducing  the  overall 
impulse. 
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We  also  note  that  Martec’s  Chinook  code,  which  makes  the  same  homogeneous  mixture 
assumption,  also  over-predicts  the  plate  loading,  Meanwhile,  the  AUTODYN  model  addresses 
the  problem  of  the  soil  becoming  a  multiphase  cloud  in  a  different  way:  the  soil  material  fails  and 
is  removed  from  the  simulation  allowing  the  gas  to  escape.  The  result  is  that  the  soil  material  will 
have  a  reduced  momentum  due  to  the  failure  mechanism.  We  note  that  failure  may  occur  due  to  a 
dropping  pressure,  but  also  cells  may  be  removed  simply  due  to  mesh  tangling.  In  this  case,  there 
is  no  physical  justification  for  the  failure  mechanism  and  the  results  are  likely  to  be  highly 
dependent  on  the  initial  mesh  configuration.  We  also  note  that  AUTODYN  has  an  option  to  allow 
the  mass  of  the  free  nodes  to  be  included  in  the  momentum  transfer  to  the  plate.  In  this  case,  a 
disembodied  node  created  through  material  failure  travels  on  a  linear  trajectory  until  it  strikes  a 
surface  and  is  not  treated  as  a  multiphase  particle  with  associated  drag  functions.  We  could  not 
verify  whether  this  mode  was  employed  in  the  AUTODYN  simulations  that  are  included  in  our 
comparisons. 

Furthermore,  there  is  some  uncertainty  in  how  to  define  the  material  properties  of  the  solids 
that  are  in  the  sandy  clay  soil.  The  prairie  soil  composition  used  in  the  MIP  experiment  consisted 
of  a  composition  of  40-60%  sand  and  10%-20%  clay.  The  exact  elastic  modulus  of  the  clay 
particles  are  not  known  to  the  investigators  at  this  time,  and  may  have  some  influence  on  the 
impulse  loading  on  the  plate.  In  the  Martec  report  [Martin  and  Link,  2003],  the  Chinook  code 
was  used  to  study  this  sensitivity  where  they  reduced  the  solid  material  elastic  modulus  by  a 
factor  of  1000.  The  code  obtained  good  agreement  with  this  value.  However,  in  our  view,  it  is 
unlikely  that  the  clay  would  cause  such  a  significant  change  in  the  elastic  modulus  of  the  solids, 
and  so  the  result  would  seem  to  suggest  that  the  impulse  loading  would  be  rather  insensitive  to  a 
small  change  (e.g.  less  than  an  order  of  magnitude)  in  the  material  stiffness. 
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Figure  3.38:  Comparison  of  impulse  on  the  target  plate  for  DOB=50mm. 
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Figure  3.39:  Comparison  of  impulse  on  the  target  plate  for  DOB=100mm. 


Figure  3.40:  Time  development  of  impulse  of  landmine  explosion  impacted  onto  a  plate  for 

DOB=50mm  and  18.2%  moisture  content. 
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Figure  3.41:  Density  contour  of  landmine  explosion  impacted  onto  a  plate  for  DOB=50mm 

and  18.2%  moisture  content  at  t=900ps. 


Figure  3.42:  Density  contour  of  landmine  explosion  impacted  onto  a  plate  for  DOB=50mm 

and  18.2%  moisture  content  at  t=1700ps. 


3.7  Observations  for  further  investigation 

From  the  small-scale  dry  sand  validation  of  the  crater  and  ejecta  cloud  evolution,  we  have 
some  evidence  that  material  strength  may  play  an  important  role  in  the  later  times  of  the  soil-blast 
progression.  The  algorithms  that  we  currently  employ  have  been  developed  primarily  for  accurate 
gas-dynamics  simulation;  therefore  including  strength  presents  additional  challenges. 
Specifically,  in  the  current  model,  we  track  material  velocity  rather  than  material  deformation, 
while  strain  is  not  directly  available.  However,  development  of  strength  models  in  the  context  of 
a  Godunov  type  solver  similar  to  the  algorithm  employed  here  have  been  developed  [Miller  and 
Colella,  2001]  and  could  be  adapted  to  our  model  with  modest  effort. 
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Probably  the  most  critical  feature  that  is  lacking  in  our  simulation  is  a  more  correct  modeling 
of  the  transition  to  multiphase  flow  of  soil  during  the  explosive  interaction.  In  the  next  quarter  we 
plan  to  implement  and  evaluate  the  effects  of  coupling  our  present  model  with  a  Lagrangian 
particle  tracking  model.  Transition  to  multiphase  would  occur  based  on  criterions  similar  to 
failure,  e.g.  if  soil  material  pressure  (or  density)  dropped  below  a  given  threshold,  we  can  assume 
that  the  soil  decohesion  results  in  a  dense  multiphase  cloud.  We  expect  that  inclusion  of  these 
effects  will  improve  our  model  performance  on  the  MIP  validation  case. 


70 


UNCLASSIFIED 


Section  4:  Simulation  of  blast  interactions  with  debris 

In  this  section,  we  investigate  the  effects  of  blast  interactions  with  debris.  One  of  the  primary 
objectives  of  this  effort  was  to  characterize  the  effects  of  these  interactions  in  order  to  develop 
appropriate  phenomenological  models  to  couple  a  multiphase  model  of  soil  particulates  with  the 
blast  event.  For  this  investigation,  we  limited  our  study  to  the  blast  effects  on  spherical  quartz 
particles  in  different  configurations  in  order  to  determine  the  effects  of  particle  configuration  on 
the  loads  that  are  communicated  to  the  particles.  We  performed  these  simulations  using  an  the 
implicit  time  integrator  in  Loci/CHEM  which  supports  over-set  meshes  and  independent 
movement  of  bodies  subject  to  a  six  degree  of  freedom  (6DOF)  rigid  body  motion.  We  first 
performed  a  validation  study  for  an  isolated  sphere  subject  to  an  impinging  shock  for  which 
experimental  data  is  available.  We  used  this  data  to  validate  the  over-set  methodology  and  also  to 
identify  the  mesh  spacing  that  will  be  required  to  obtain  a  reasonable  resolution  of  particle  loads. 
We  employed  simulations  of  an  isolated  sphere  subjected  to  a  strong  shock  at  various  Reynolds 
numbers  to  determine  the  basic  requirements  of  the  simulation.  We  then  considered  6DOF 
models  for  various  particle  configurations  including  tandem  particle  arrangements  with  different 
particle  spacings  as  well  as  staggered  configurations.  Finally,  we  performed  a  6DOF  simulation 
on  a  collection  of  particles  that  interact  with  a  spherical  blast  wave  in  place  of  the  planar  shock. 

4.1  Shock-sphere  validation  case 

Our  time  dependent,  over-set  mesh  test  case  is  the  comparison  of  experimental  drag  forces 
computed  on  a  sphere  that  is  subjected  to  a  passing  shock  wave.  The  experimental  setup, 
performed  at  the  Institute  of  Fluid  Science,  Tohoku  University,  Japan,  is  described  in  a  paper  by 
Tanno  et.  al.  [Tanno  et  al.,  2004],  In  this  experiment,  a  sphere  was  subjected  to  aM=1.22  shock 
wave  impingement.  The  geometry  for  this  case  consisted  of  a  300mmx300mm  square  cross 
section  shock  tube  where  an  80 mm  diameter  sphere  was  suspended  on  an  elastic  wire.  The  driver 
side  of  the  shock  tube  consisted  of  nitrogen  gas  at  a  pressure  of  P=240kPa  at  room  temperature, 
while  the  driven  gas  was  air  at  atmospheric  conditions.  The  post  shock  conditions  for  this  case 
were  a  pressure  of  P=\59.0SkPa,  a  static  temperature  T=334K,  and  a  velocity  of  1 1 4.6mA.  The 
force  exerted  on  the  sphere  as  the  shock  passed  was  measured  through  accelerometers  embedded 
in  the  sphere.  For  the  duration  of  the  experiment,  the  sphere  movement  was  slight  and  one  can 
reasonably  approximate  the  sphere  as  a  rigid  body  for  computing  the  drag  force  for  comparison 
with  the  experimental  data. 

We  wish  to  validate  the  over-set  mesh  methodology  that  we  will  employ  for  our  many-body 
problem  to  facilitate  including  debris  motion.  In  this  case,  we  are  validating  the  utility  of  using 
this  over-set  methodology  for  accurately  capturing  shock-body  interactions.  In  addition,  we  are 
attempting  to  get  an  estimate  of  the  mesh  resolution  required  to  obtain  a  reasonably  accurate 
solution.  The  more  coarse  the  mesh  we  can  use,  potentially  the  more  independent  bodies  we  can 
consider  in  our  debris  simulation. 

In  this  case  we  constructed  the  over-set  mesh  using  a  Cartesian  mesh  matching  the  square 
cross  section  of  the  shock  tube  and  used  an  unstructured  mesh  to  capture  the  geometry  of  the 
sphere.  For  the  sphere  mesh,  we  generated  an  unstructured  mesh  that  extended  to  an  interface 
sphere  with  a  radius  of  120mm.  A  Cartesian  background  mesh  was  generated  for  the 
1000mmx300mmx300mm  shock  tube  section.  Mesh  sizes  were  generated  with  a  mesh  spacing  of 
10mm,  5mm,  and  2.5mm  with  the  combined  mesh  sizes  being  0.264 M,  1.64 M,  and  11.5 M  cells 
respectively.  For  the  all  of  these  cases,  a  time-step  of  l|xs  and  5  Newton  iterations  with  second 
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order  time  differencing  were  employed.  The  Reynolds  number  for  this  case  was  7.6x1 06. 
Menter’s  Shear  Stress  Transport  (SST)  turbulence  model  was  utilized  with  solid  viscous  walls 
modeled  using  a  wall  law  boundary  condition  with  a  viscous  wall  boundaries  meshed  with  a 
normal  spacing  giving  y+  in  the  range  of  50-100. 

In  Figure  4.1,  the  details  of  the  mesh  topology  are  shown  for  the  5 mm  case  along  with  the 
density  contours  at  a  time  200ps  after  the  shock  hits  the  sphere.  The  cutting  plane  shows  a  cross 
section  of  the  mesh  illustrating  the  Cartesian  mesh  of  the  shock  tube  and  the  unstructured  mesh 
that  surrounds  the  sphere.  Note  that  the  unstructured  sphere  mesh  is  automatically  cut  by  the 
over-set  mesh  facility  due  to  proximity  with  the  walls  of  the  shock  tube.  Interpolation  is  used  to 
couple  the  solution  at  the  interface  of  the  Cartesian  and  unstructured  meshes. 

The  comparisons  of  the  simulation  results  with  the  experimental  data  are  shown  in  Figure  4.2. 
All  grids  produced  reasonable  comparisons  to  the  experimental  results,  with  the  5 mm  and  2.5mm 
grids  producing  nearly  identical  results.  With  the  exception  of  a  small  peak  in  the  experimental 
data  at  around  50ps,  the  simulations  on  all  mesh  resolutions  produced  excellent  agreement  with 
the  experimental  data.  In  addition,  both  experiment  and  simulations  reproduced  the  curious  effect 
of  negative  drag  (occurring  around  400j.lv)  that  is  caused  by  a  coalescence  of  shocks  on  the  back 
side  of  the  sphere.  This  test  case  shows  that  the  code  can  accurately  capture  shocks  passing 
though  these  interpolated  regions.  We  also  note  that  the  10 mm  mesh  spacing,  although  very 
coarse,  provided  remarkably  good  agreement  with  experimental  data.  This  suggests  that  we  may 
be  able  to  sacrifice  mesh  resolution  without  sacrificing  the  accuracy  of  our  prediction  in  order  to 
simulate  a  large  number  of  bodies. 


Figure  4.1:  Density  contours  on  cutting  plane  showing  over-set  mesh. 
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Drag  Force  on  a  Shocked  Sphere 


4.2  Study  of  shock  wave  impinging  on  a  single  sphere 

Before  we  performed  the  6DOF  simulations  of  blast-particle  interactions,  we  first  performed 
simulations  of  an  isolated  particle  of  various  sizes  with  a  range  of  impingement  shock  strengths 
in  order  to  understand  the  time  scales  of  various  phenomena.  Irregularity  of  the  particle  shape  is 
beyond  the  scope  of  the  current  study;  thus,  only  solid  particles  with  spherical  shape  were 
considered.  The  particle  was  placed  in  a  shock  tube  and  the  time  history  of  the  momentum  flux 
exerted  on  the  particle  was  computed. 

Many  studies  have  been  conducted  on  shock  wave  interactions  with  particles  at  different 
length,  time,  and  velocity  scales  [Sun  et  al.,  2004,  Liou  and  Takayama,  2004,  Tanno  et  al.,  2003, 
Saito  et  ah,  2003],  The  dynamic  drag  coefficient  of  a  sphere  due  to  shock  wave  loading  was 
investigated  numerically  and  experimentally  by  Sun  et.al  [Sun  et  ah,  2004],  In  their  study,  the 
unsteady  drag  of  a  sphere  due  to  a  shock  wave  of  Mach  number  1 .22  was  examined  over  a  range 
of  Reynolds  numbers  varying  from  102  to  105,  which  is  equivalent  to  the  diameter  of  the  sphere 
varying  from  8 \xm  to  80 mm.  It  was  found  that  unsteady  pressure  loading  due  to  shock  interaction 
with  the  sphere  is  a  dominant  effect  early  in  the  interaction,  especially  for  high  Reynolds  number 
cases.  The  effects  of  viscosity  on  the  drag  are  negligible  at  this  stage,  although  they  become  more 
important  after  the  shock  passes.  The  duration  of  the  strong,  unsteady  shock-particle  interaction 
is  approximately  the  time  for  the  shock  wave  to  travel  one  sphere  diameter.  In  our  study,  we 
intended  to  determine  reasonable  parameters  for  particles  (e.g.  particle  diameter)  and  fluid 
conditions  (e.g.  Mach  number)  so  that  we  have  a  basic  understanding  of  how  particles  will 
respond  during  the  time  when  the  shock  passes.  Analysis  based  on  the  momentum  equation  for 
the  particle  showed  that  the  speed  the  particle  attains  due  to  the  initial  shock  impingement  does 
not  depend  on  the  size  of  the  particle,  but  rather  on  the  density  of  the  particle,  the  Mach  number 
of  the  flow,  and  the  density  of  the  gas  behind  the  approaching  shock.  According  to  our  analysis, 
for  quartz,  with  a  density  of  2650 kg/m3  in  our  application,  particles  will  have  little  response 
during  the  short  duration  of  time  that  corresponds  to  the  shock  wave  passing  over  the  particle. 
Thus,  the  response  of  the  particle  to  the  initial  blast  wave  will  be  dependent  on  the  strength  of  the 
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wave  and  the  relative  densities  of  the  blast  wave  and  the  material  density  of  the  particle. 
Therefore,  unless  the  particle  is  very  close  to  the  expanding  explosive  cloud,  the  particle  response 
will  be  much  slower  than  the  approaching  shock.  It  is  important  to  note  that  this  observation 
appears  to  be  true  regardless  of  the  particle  size  provided  that  the  Reynolds  number  of  the  particle 
is  on  the  order  of  103  or  greater.  Over  much  longer  durations  of  time,  the  effects  of  particle  sizes 
will  become  more  important. 

Movement  of  the  spherical  particles  due  to  the  shock  wave  was  computed  using  6DOF,  where 
no  model  was  used  for  the  coupling  between  flow  and  the  particles.  The  above  speculations  on 
particle  response  to  the  flow  will  be  tested.  As  preparation  for  the  6DOF  simulation,  the  drag  on 
stationary  particles  with  different  diameters  under  different  shock  speeds  were  simulated  and 
compared  with  analysis. 

Two  Mach  numbers,  2.5  and  5.0,  were  chosen  for  the  impinging  shock  strengths  and  three 
different  particle  sizes  were  considered  for  each  Mach  number.  Three  representative  particle 
diameters  in  soil  are  selected  to  be  0.05 mm,  0.5 mm  and  5 mm,  which  cover  the  range  of  sizes  from 
a  typical  sand  grain  to  small  gravel.  The  Reynolds  number  ranges  from  3xl03  to  3xl05  for  Mach 
2.5,  and  8xl03  to  8xl05  for  Mach  5.  The  definition  of  Reynolds  number  is 

„  P«U«d 

Re=——,  (4.1) 

Ha 


where  d  is  the  diameter  of  the  sphere.  This  definition  of  the  Reynolds  number  has  been  generally 
adopted  to  document  the  drag  coefficients  for  a  sphere,  p  ,  U ,  pQ  are  characteristic  density, 

velocity  and  viscosity  of  the  flow,  and  are  taken  as  those  behind  the  shock  wave. 

The  results  from  our  simulations  are  illustrated  in  Figures  4.3  and  4.4.  These  two  plots  are  the 
dimensionless  time  history  of  the  non-dimensional  drag  force  (CJ  for  two  different  Mach 

numbers.  The  physical  time  is  non-dimensionalized  using  r/U  ,  where  r  is  the  radius  of  the  sphere 
and  U  is  the  post-shock  fluid  velocity.  The  non-dimensional  drag  force  is  the  drag  coefficient, 
which  is  defined  as 


c=  y 

d  Po^r2’ 


(4.2) 


where /is  the  drag  force  on  the  sphere.  These  simulations  show  that  the  particle  size  has  minimal 
effect  on  the  particle  response  to  the  shock  wave  in  these  size  categories.  Note  that  these  curves 
overlay  one  another  because  of  the  non-dimensionalization  applied.  Thus,  a  particle  is  exposed  to 
the  strong  shock  loads  for  a  shorter  period  of  time  for  smaller  size  particles,  while  the  forces  are 
also  reduced  by  the  relative  difference  in  the  projected  surface  area  of  the  particles.  As  a  result, 
although  a  larger  particle  will  receive  an  increased  force,  this  force  will  be  compensated  by  its 
increased  mass  such  that  the  resulting  acceleration  the  particle  experiences  will  be  similar 
regardless  of  particle  size.  This  result  suggests  that  for  particles  that  have  mass  similar  to  soil, 
there  will  be  very  little  motion  of  the  particles  during  the  immediate  interaction  with  the  blast 
wave,  regardless  of  the  particle  size. 
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Figure  4.3:  Time  development  of  the  drag  coefficient  on  spheres  with  impinging  shock 

speed  of  M=2.5. 


Figure  4.4:  Time  development  of  the  drag  coefficient  on  spheres  with  impinging  shock 

speed  of  M=5. 
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4.3  Study  of  a  normal  shock  wave  impinging  on  multiple  spheres 

In  this  section,  we  will  discuss  the  interaction  of  the  flow  field  with  multiple  spheres.  First, 
we  study  the  flow  of  four  spheres  in  a  tandem  configuration,  separated  by  an  equal  distance,  as 
illustrated  in  Figure  4.5.  Three  sphere  diameters  are  taken  to  be  5 mm,  0.5 mm  and  0.05  mm.  For 
each  sphere  diameter,  simulations  are  performed  for  three  equal  inter-sphere  spacings  of  s=\r, 
s=2r  and  s=6r,  as  illustrated  in  Figure  4.5.  A  shock  wave  moving  at  Mach  5.0  relative  to  the 
ambient  air  impinges  on  the  spheres.  Based  on  the  previous  validation  exercise,  a  mesh  size  of 
0.0625  times  the  sphere  diameter  is  adopted  in  the  simulations  for  each  sphere  size. 


Figure  4.5:  Schematic  diagram  of  tandem  particle  configuration. 

One  of  the  important  parameters  for  the  study  is  the  particle  drag  coefficient  as  defined  in 
equation  (4.2).  In  these  cases,  pQ  and  U  are  the  characteristic  density  and  velocity,  respectively, 

and  are  taken  to  be  the  ambient  density  before  the  shock  arrives  and  the  velocity  of  the  shock 
relative  to  the  ambient  velocity.  According  to  equation  (4.1),  the  Reynolds  numbers  are  6.4x1 05, 
6.4xl04  and  6.4xl03  for  5 mm,  0.5 mm  and  .0.5 mm  diameter  spheres,  respectively.  An  example  of 
the  simulated  pressure  contours  for  a  particle  size  of  5 mm  and  a  shock  speed  of  M=  5  with  s=2r  is 
shown  in  Figure  4.6. 


Figure  4.6:  Simulated  pressure  contours  for  M=5,  r=2.5mm,  and  s=2r. 

The  results  of  time  evolution  of  the  drag  coefficient  on  the  particles  are  shown  in  Figures  4.7- 
4.15  for  three  different  sphere  diameters  at  different  spacings.  Since  there  is  a  wide  range  of 

r 

length  scales  for  the  spheres,  the  physical  time  is  normalized  by  jt,  and  non-dimensional  time  is 

uo 

used  for  the  images.  The  figures  show  that  for  the  same  spacing  between  spheres,  the  drag 
coefficients  are  almost  the  same  for  the  three  different  particle  sizes.  This  indicates  that  pressure 
drag  (inviscid  force)  is  the  main  contributor  to  the  total  drag,  and  the  viscous  drag  is  not 
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significant  for  the  particle  Reynolds  numbers  studied  in  the  present  cases.  It  is  also  evident  that 
the  spheres  experience  a  steep  drag  increase  due  to  shock  impingement  and  that  the  drag  drops 
precipitously  as  the  shock  wave  moves  to  the  rear  side  of  the  sphere.  The  maximum  drag  on 
spheres  2,  3  and  4  is  larger  than  that  on  the  first  sphere  because  the  flow  accelerates  after  passing 
the  leading  sphere  and  is  then  focused  on  the  following  sphere.  The  shock  weakens  as  it  travels 
downstream  with  the  presence  of  the  spheres  because  of  the  effects  of  shock  reflections,  which 
may  be  considered  blockage  effects.  The  blockage  effect  is  more  obvious  when  the  spheres  are 
closer  together;  thus,  the  maximum  drag  exerted  on  spheres  3  and  4  decreases  for  1  r  and  2 r 
spacing,  as  shown  in  the  figures.  However,  when  the  distance  between  spheres  is  larger,  the 
blockage  effect  becomes  less  significant,  as  does  the  focusing  effect.  Therefore,  the  spheres 
behave  more  like  isolated  individual  particles  without  much  interaction,  as  depicted  in  Figures 
4.9,  4.12,  and  4.15.  Another  important  observation  made  from  those  figures  is  that  the  drag  values 
at  steady  state  for  the  first  and  last  spheres  are  higher  than  those  for  the  middle  spheres.  There  are 
two  phenomena  that  are  responsible  for  this  effect:  1)  the  leading  sphere  experiences  an  increased 
pressure  on  its  trailing  side  due  to  a  supporting  effect  of  the  following  sphere,  while  the  trailing 
sphere  experiences  a  much  lower  pressure  on  the  trailing  side  and  2)  the  leading  sphere  sees  an 
increased  stagnation  pressure,  while  the  trailing  spheres  experience  a  shielding  effect  from  the 
spheres  ahead  of  it.  So  the  middle  spheres  experience  both  a  shielding  effect  from  the  particle 
ahead  of  it  and  a  supporting  effect  from  the  particle  that  trails  it.  The  leading  particle  does  not 
benefit  from  a  shielding  effect,  while  the  trailing  particle  does  not  benefit  from  a  supporting 
effect.  Therefore,  we  expect  more  loading  will  be  experienced  at  the  leading  and  trailing  edges  of 
the  particulate  cloud.  As  seen  in  the  figures,  the  difference  between  steady  state  drag  on  spheres 
at  different  locations  appears  to  be  smaller  as  the  spacing  between  the  spheres  becomes  larger 
since  the  shielding  and  supporting  effects  related  to  the  particle-particle  interactions  tend  to 
decrease  with  increasing  particle  spacing. 


Figure  4.7:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
5 mm  diameter  and  v=  1  r. 
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Figure  4.8:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
S/;;/;/  diameter  and  s=2r. 


Figure  4.9:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
5/77/7/  diameter  and  s=6r. 
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Figure  4.10:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.5 mm  diameter  and  s=  1  r. 


Figure  4.11:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.5mm  diameter  and  s=2r. 
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Figure  4.12:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.5 mm  diameter  and  s=6r. 


Figure  4.13:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.05mm  diameter  and  s=  1  r. 
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Figure  4.14:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.05mm  diameter  and  s=2r. 


Figure  4.15:  Drag  coefficient  of  shock  wave  loading  on  spheres  with 
0.05mm  diameter  and  s=6r. 


Figures  4.16-4.19  compare  the  drag  coefficient  for  the  three  different  diameters  of  the  spheres 
at  the  largest  particle  spacing  ( s=6r ).  As  mentioned  above,  the  particle  drag  coefficients  for  each 
different  size  of  sphere  are  almost  the  same  for  the  particle  Reynolds  number  in  our  present 
studies,  which  is  again  demonstrated  in  this  plot.  Although  the  viscous  drag  in  not  significant,  it 
is  still  observed  that  the  viscous  effect  increases  with  the  decrease  of  particle  Reynolds  number 
(smaller  particle  size). 
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Figure  4.16:  Comparison  of  drag  coefficients  on  leading  sphere  with 
different  diameters  for  s=6r. 


Figure  4.17:  Comparison  of  drag  coefficients  on  second  sphere  with 
different  diameters  for  s=6r. 
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Figure  4.18:  Comparison  of  drag  coefficients  on  third  sphere  with 
different  diameters  for  s=6r. 


Figure  4.19:  Comparison  of  drag  coefficients  on  trailing  sphere  with 
different  diameters  for  s=6r. 

Figures  4.20-4.22  show  the  computed  particle  velocities  for  the  5mm  diameter  tandem 
spheres.  The  velocities  in  these  plots  were  shifted  so  that  they  all  start  with  a  zero  time  at  the  time 
of  arrival.  If  the  inter-particle  interactions  are  negligible,  we  would  expect  these  curves  to  be 
identical.  We  note  for  a  sphere  spacing  of  1  sphere  radius,  we  have  significant  differences 
between  the  particles  with  the  first  and  last  sphere  gaining  the  most  momentum  from  the  flow. 
The  middle  particles  acquired  significantly  less  momentum  from  the  shock  impingement.  The 
middle  particles  have  both  a  shielding  effect  from  the  wake  of  the  sphere  ahead  of  it,  while  they 
also  have  a  supporting  effect  from  their  own  wake  interacting  with  the  sphere  that  trails  them. 
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Both  of  these  effects  have  the  result  of  reducing  the  drag  on  the  particle.  Since  the  leading 
particle  does  not  have  a  shielding  effect,  and  the  trailing  particle  does  not  have  a  supporting 
effect,  these  two  particles  experience  a  larger  drag  than  the  middle  particles.  We  see  that  as  we 
increase  the  sphere  spacing  to  1  sphere  diameter,  the  degree  of  difference  between  the  leading, 
middle,  and  trailing  particles  is  much  less,  as  is  observed  in  Figure  4.21.  For  the  case  where  the 
sphere  spacing  is  3  diameters,  as  shown  in  Figure  4.22,  we  see  that  the  spheres  experience  a 
decreasing  momentum  transfer  as  we  proceed  down  the  line  of  particles.  However,  the  leading 
and  trailing  particles  do  not  show  a  particular  departure  from  the  trend.  We  observe  that  this 
effect  could  simply  be  explained  by  the  loss  of  fluid  momentum  due  to  the  drag  of  the  particle 
ahead  of  it,  and  the  shielding  and  supporting  effects  appear  to  no  longer  be  playing  a  role  in  the 
fluid-particle  coupling. 

The  variation  of  sphere  velocities  with  the  spacing  between  them  is  plotted  in  Figures  4.23- 
4.26.  It  is  expected  that  the  velocity  increases  with  increasing  sphere  spacing  since  the  damping 
(including  the  shielding  and  supporting  effect)  from  other  spheres  weakens  with  increasing 
distance  between  the  spheres.  The  velocity  of  leading  and  middle  spheres  (Figures  4.23,  4.24,  and 
4.25)  supports  this  hypothesis,  while  the  reverse  trend  is  observed  for  the  trailing  sphere.  The 
anomaly  may  come  from  the  complex  behavior  in  the  wake  region.  It  is  seen  that  the  velocities  of 
the  two  middle  spheres  are  more  sensitive  to  the  spacing  since  they  feel  both  the  shielding  and 
supporting  effects  from  the  spheres  in  front  and  behind,  while  the  first  and  last  spheres  do  not 
experience  as  much  interaction  from  other  spheres.  We  also  note  that  the  interactions  appear  to 
be  complex  and  nonlinear.  For  example,  in  Figure  4.23,  the  lowest  particle  velocity  is  obtained 
for  the  intermediate  spacing.  From  the  graph  it  appears  that  this  is  due  to  a  transient  interaction 
between  particles  that  occurs  at  about  10  microseconds.  Nonetheless,  the  general  slope  after  this 
transient  event  still  supports  our  hypothesis  regarding  the  supporting  effect  on  the  leading  sphere. 


Figure  4.20:  Comparison  of  velocity  for  different  spheres  ( d=5mm ,  s=lr). 
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Figure  4.21:  Comparison  of  velocity  for  different  spheres  ( d=5mm ,  s=2r). 


Figure  4.22:  Comparison  of  velocity  for  different  spheres  ( d=5mm ,  s=6r). 
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time  (s) 

Figure  4.23:  Comparison  of  velocity  of  leading  sphere  ( d=5mm ) 
for  different  particle  spacings. 


time  (s) 

Figure  4.24:  Comparison  of  velocity  of  second  sphere  (d=5mm) 
for  different  particle  spacings. 
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Figure  4.25:  Comparison  of  velocity  of  third  sphere  ( d=5mm ) 
for  different  particle  spacings. 


Figure  4.26:  Comparison  of  velocity  of  trailing  sphere  ( d=5mm ) 
for  different  particle  spacings. 

The  predicted  movement  of  the  spheres  is  illustrated  in  Figures  4.27-4.29.  The  four  spheres 
travel  a  small  distance  and  do  not  have  the  chance  to  collide  during  the  time  frame  simulated  in 
the  present  study.  Figures  4.30-4.33  compare  the  movement  of  spheres  at  different  sphere 
spacings.  The  distance  the  spheres  move  is  in  agreement  with  the  velocity  profile  depicted  in 
Figures  4.23-4.26,  but  the  spheres  with  larger  spacing  start  moving  at  a  later  time,  as  expected. 

To  reduce  the  computational  expense,  all  the  above  results  were  obtained  in  the  one-way 
coupled  frame  in  which  the  spheres  are  stationary  and  the  time-integrated  drag  force  is  employed 
to  compute  the  velocities  and  positions  as  a  post-processing  procedure.  However,  one  simulation 
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was  performed  in  the  two-way  coupled  frame  where  the  spheres  move  in  the  flow  and  exert  an 
influence  on  the  flow  field.  The  sphere  moves  with  six  degrees  of  freedom  (6DOF)  under  the 
integrated  force  from  the  flow  on  the  sphere  surface  and  sphere  body  force.  Figure  4.34  compares 
the  results  between  the  one-way  and  two-way  coupled  simulations.  The  two  types  of  simulations 
provide  identical  results  for  the  current  problem. 


Figure  4.27:  Position  of  spheres  ( d=5mm )  along  shock  tube  axis  with  v=  1  r. 


time  (s) 

Figure  4.28:  Position  of  spheres  (d=5mm)  along  shock  tube  axis  with  s=2r. 
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Figure  4.29:  Position  of  spheres  ( d=5mm )  along  shock  tube  axis  with  s=6r. 


Figure  4.30:  Position  of  leading  sphere  ( d=5mm )  along  shock  tube 
axis  for  different  spacings. 
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Figure  4.31:  Position  of  second  sphere  ( d=5mm )  along  shock  tube 
axis  for  different  spacings. 


Figure  4.32:  Position  of  third  sphere  ( d=5mm )  along  shock  tube 
axis  for  different  spacings. 
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Figure  4.33:  Position  of  trailing  sphere  (d=5mm)  along  shock  tube 
axis  for  different  spacings. 
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Figure  4.34:  Comparison  of  positions  for  stationary  and  moving  spheres. 
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4.4  Study  of  a  normal  shock  wave  impinging  on  multiple  staggered 
spheres 

Next,  a  staggered  sphere  configuration  was  considered.  This  configuration  is  an  arrangement 
of  six  spheres,  with  one  leading  sphere,  four  middle  spheres,  and  one  trailing  sphere.  A  diagram 
of  the  layout  of  the  four  spheres  that  are  aligned  with  the  z=0  plane  is  shown  in  Figure  4.35. 
Sphere  #3  and  sphere  #5  are  not  shown  and  will  appear  at  positive  and  negative  z  coordinates, 
respectively.  The  lines  at  the  top  and  bottom  of  this  schematic  are  the  location  of  periodic 
boundaries. 


Figure  4.35:  Schematic  diagram  of  sphere  layout  for  staggered  case  showing  spheres 

located  on  z= 0. 

Figure  4.36  shows  an  instantaneous  snapshot  of  pressure  contours  just  as  the  Mach  5.0  shock 
wave  is  passing  over  the  trailing  sphere  in  the  staggered  configuration.  Figure  4.37  shows  the 
time  history  of  the  drag  coefficient  on  the  six  spheres.  The  results  on  the  four  spheres  in  the 
middle  plane  are  almost  identical  due  to  the  symmetry  of  the  particle  layout.  The  maximum  drag 
on  the  four  middle  spheres  is  almost  equivalent  to  the  drag  on  the  first  one  because  there  are  no 
spheres  in  front  directly  blocking  the  flow,  which  provides  a  shielding  effect.  In  addition,  we 
don’t  observe  a  significant  focusing  effect  due  to  the  presence  of  the  leading  sphere.  A  sharp 
increase  of  maximum  drag  on  the  last  sphere  (sphere  6)  is  caused  by  the  fact  that  the  flow 
accelerates  when  passing  the  middle  four  spheres,  which  focuses  the  flow  down  the  central  core 
directly  onto  the  trailing  sphere.  Contrary  to  the  tandem  case,  the  steady  state  drag  on  the  middle 
four  spheres  is  higher  since  there  are  no  shielding  and  supporting  effects  presented  to  those 
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spheres.  The  lower  steady  state  drag  force  on  the  last  sphere  results  from  the  shielding  effect 
from  the  first  sphere. 

The  velocities  and  positions  of  the  six  spheres  in  the  staggered  configuration  are  examined  in 
Figures  4.38  and  4.39.  As  expected,  the  middle  four  spheres  move  faster  when  compared  to  the 
first  and  last  spheres  since  there  is  no  direct  impact  from  other  spheres.  Like  the  tandem  case,  the 
spheres  still  do  not  travel  fast  enough  to  collide  with  each  other. 


Figure  4.36:  Instantaneous  pressure  contours  for  staggered  configuration  ( d=5mm ,  s=lr). 


Figure  4.37:  Drag  coefficient  on  spheres  for  staggered  configuration  ( d=5mm ,  s=lr). 
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Figure  4.38:  Velocity  of  spheres  ( d=5mm ,  s=lr)  for  staggered  configuration. 


Figure  4.39:  Distance  of  travel  for  spheres  ( d=5mm ,  s=lr)  for  staggered  configuration. 

In  order  to  demonstrate  how  the  layout  of  particles  can  affect  the  particle  dynamics,  a 
comparison  of  velocity  and  position  between  the  staggered  and  tandem  configurations  is  provided 
in  Figure  4.41.  Only  the  middle  spheres,  having  the  same  distance  (in  the  axial  shock  tube 
direction)  relative  to  the  first  sphere,  are  selected  for  this  comparison.  It  is  obvious  that  the 
particle  configuration  plays  a  significant  role  in  the  dynamics  of  spheres.  The  shielding  and 
supporting  effects  on  the  middle  sphere  in  the  tandem  layout  reduces  the  movement  of  the  sphere 
significantly. 
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Figure  4.40:  Comparison  of  velocity  of  spheres  ( d=5mm ,  s=lr)  between 
staggered  and  tandem  configurations. 


Figure  4.41:  Comparison  of  distance  traveled  for  spheres  ( d=5mm ,  s=lr)  between 

staggered  and  tandem  configurations. 

From  the  numerical  studies  of  the  shock  interaction  with  micro-spheres  for  a  particle  Reynolds 
number  ranging  from  6xl03  to  6xl05,  we  may  draw  the  following  conclusions.  First,  the  size  of 
the  particles  in  the  regime  of  the  current  application  does  not  affect  the  dynamics  of  the  particles. 
Second,  the  spacing  between  multiple  particles  and  the  layout  of  the  particles  play  an  important 
role  in  the  movement  of  the  particles.  The  shielding  effect  from  upstream  and  the  supporting 
effect  from  downstream  particles  can  slow  down  the  particle,  and  these  effects  are  most 
pronounced  in  the  tandem  configuration.  However,  the  focusing  effect  from  the  acceleration  of 
the  flow  past  the  sphere  could  impose  a  boosting  force  on  the  sphere  downstream.  Finally,  the 
unsteady  force  from  the  shock  impingement  does  not  have  a  significant  influence  on  particle 
dynamics.  Instead,  the  steady  state  force  plays  a  much  more  significant  role. 
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4.5  Study  of  a  blast  wave  impinging  on  multiple  particles 

After  investigating  the  behavior  of  particles  in  a  shock  tube,  we  conducted  a  study  on  a 
notional  blast  wave  interaction  with  particles.  To  generate  the  blast  wave,  we  performed  an 
axisymmetric  simulation  of  a  surface-laid  hemispherical  charge  of  TNT  with  a  radius  of  144mm 
and  a  mass  of  10.19 kg.  We  used  an  ideal  gas  model  with  high-pressure  conditions  within  a  sphere 
to  model  the  detonated  TNT.  The  two-dimensional  simulation  was  run  until  the  blast  wave  had 
expanded  to  about  7=1. 8m  from  the  center  of  the  TNT  charge.  At  this  point,  we  interpolated  the 
axisymmetric  solution  onto  a  three-dimensional  grid  that  consisted  of  a  solid  angle  slice  of  a 
spherical  domain.  We  used  reflecting  conditions  on  the  boundaries  of  this  solid  angle.  The  three 
dimensional  grid  in  which  we  embedded  the  particles  is  shown  in  Figure  4.42.  Before  embedding 
the  particles,  we  utilized  adaptive  mesh  refinement  in  the  regions  where  we  would  be  resolving 
the  blast-particle  interaction.  Specifically,  we  refined  two  levels  over  the  general  region  where 
the  blast  wave  would  be  present  for  the  simulation  and  another  two  levels  in  the  immediate 
region  where  the  particles  were  placed,  as  illustrated  in  Figure  4.43.  Ten  particles  were  arranged 
in  a  staggered  formation  configured  in  two  planes.  Each  plane  contains  a  configuration  similar  to 
a  “five”  die  face  with  the  subsequent  plane  rotated  by  45  degrees.  These  ten  particles  were  given 
a  diameter  of  1cm  and  density  of  2650/cg/m3  and  were  located  just  before  the  front  of  the 
interpolated  blast  wave.  These  ten  particles  are  equally  divided  into  two  planes  located  at 
7=1. 82m  and  7=1. 84m  perpendicular  to  the  approaching  blast  wave,  as  illustrated  in  Figure  4.44. 

Figures  4.45-4.49  show  the  evolution  of  the  blast  wave  impinging  on  the  particle 
configuration.  In  this  figure,  symmetry  illustrates  how  the  simulation  represents  a  larger  particle 
cloud.  From  this  simulation,  a  complex  set  of  shock  structures  emerges  that  is  similar  to  the 
normal  shock  simulations  described  earlier.  In  Figure  4.49,  there  is  considerable  variation  in 
particle  velocities,  as  illustrated  by  the  particle  colors.  The  source  of  these  differences  in  particle 
velocity  can  be  seen  in  the  time  history  of  the  particle  drag  forces,  as  shown  in  Figures  4.50  and 
4.51.  Compared  with  particles  in  the  shock  wave  presented  in  the  last  section,  the  drag  force  on 
spheres  in  blast  waves  exhibit  similar  general  behavior  of  peaking  and  falling.  However,  the  drag 
force  at  steady  state  after  the  blast  wave  passes  decreases  rapidly,  which  is  a  result  of  the  blast 
wave’s  expanding  surface  area,  causing  a  rapidly  decreasing  pressure  profile. 


96 


UNCLASSIFIED 


Figure  4.42:  Solid  angle  partition  of  sphere  used  for  blast-particle  simulation. 


Figure  4.43:  Refined  solid  angle  mesh  in  particle  region. 
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Figure  4.44:  Configuration  of  ten  particles  used  in  blast-particle  simulation. 


Figure  4.45:  Time  evolution  of  blast-particle  interaction  at  50ps. 
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Figure  4.46:  Time  evolution  of  blast-particle  interaction  at  150ps. 


Figure  4.47:  Time  evolution  of  blast-particle  interaction  at  250ps. 
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Figure  4.48:  Time  evolution  of  blast-particle  interaction  at  350ps. 


Figure  4.49:  Time  evolution  of  blast-particle  interaction  at  450ps. 
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Figure  4.50:  Time  development  of  drag  on  leading  plane  spheres  in  blast  wave. 


Figure  4.51:  Time  development  of  drag  on  trailing  plane  spheres  in  blast  wave. 
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Section  5:  Simulation  of  coupled  fluid-structure 
interactions 

A  coupled  air-blast-wave  (shock)- structure  interaction  analysis  capability  utilizing 
Loci/CHEM  and  LS-DYNA  has  been  developed  and  demonstrated  via  the  fluid- structure 
interaction  of  steel  plates  for  various  TNT  charge  weights.  Additionally,  loosely-coupled  and 
one-way  coupled  strategies  were  examined  using  a  two-dimension  blast  wave  interacting  with  an 
elastic  circular-arc  “bump”  geometry. 

5.1  Two-dimensional  FSI  simulations 

The  software  coupling  strategies  under  examination  included  a  loosely  coupled  (“two-way”) 
methodology  and  a  decoupled  (aka  “one  way”)  methodology.  The  loosely  coupled  approach 
solved  each  disciplinary  set  of  equations  within  its  own  domain  and  interfaced  the  solutions  at 
common  boundaries.  For  this  coupling  strategy,  even  when  one  system  had  a  linear  equation  set 
(e.g.,  linear  elastic  structural  analysis),  the  problem  remained  nonlinear.  For  example,  in  fluid- 
structure  interaction  (FSI)  with  nonlinear  CFD  and  linear  computational  structural  dynamics 
(CSD),  the  external  force  vector  is  an  implicit  function  of  the  displacements.  Thus,  during  each 
time  step,  dynamic  equilibrium  between  the  systems  must  be  achieved  via  sub-iterations. 
However,  situations  arise  when  a  two-way  coupling  of  the  disciplinary  physics  is  deemed 
unnecessary  and  neglecting  the  interaction  is  considered  to  be  a  valid  approximation.  This 
situation  arises  when  the  deformation  of  the  structure  has  little  impact  on  the  developing  blast 
(that  is,  when  the  load  duration  is  much  smaller  than  the  response  time  of  the  structure)  and, 
hence,  a  one-way  coupling  is  considered  sufficient. 

Regardless  of  coupling  strategy,  the  FSI  interdisciplinary  transfer  algorithm  must  be  utilized 
for  communication  between  the  disciplinary  domains,  which  typically  do  not  have  the  same 
discretizations  at  the  common  interface.  Accuracy  of  these  transfer  algorithms  may  be  evaluated 
by  comparing  the  transferred  disciplinary  data  to  those  of  the  matched-mesh  solution  (that  is, 
same  discretization  along  the  interface  for  both  disciplines).  To  this  end,  a  2D  “testbed”  tool  has 
been  developed,  henceforth  referred  to  as  2DFSI.  An  effort  is  underway  to  assess  the  accuracy  of 
the  in-house  non-overlapping  point-to-point  interdisciplinary  data  transfer  scheme  and  to 
investigate  the  disciplinary  coupling  strategies  (e.g.,  loose  “two-way”  coupling,  and  decoupled  or 
“one-way”  coupling).  A  brief  description  of  the  testbed  software  and  preliminary  results  will  be 
presented  and  discussed  within  this  report. 

Currently,  the  computer  code  used  to  obtain  the  flow  solution  in  2DFSI  is  called  NS2D.  As  the 
name  implies,  NS2D  is  an  unsteady  implicit  2-D  Navier-Stokes  flow  solver.  The  code  can 
operate  in  either  a  viscous  or  inviscid  mode.  The  convective  flux  vectors  are  discretized  using  a 
high-resolution  flux-difference-split  TVD  scheme,  and  the  diffusive  flux  vectors  are  discretized 
using  central  difference  operators.  The  same  software  was  previously  coupled  to  a  heat  energy 
conduction  code  in  an  aero-thermal  design  endeavor  [Janus  and  Newman,  2000].  The  motivation 
and  justification  for  utilizing  this  software  at  this  stage  (rather  than  Loci/CHEM)  are  that  the  code 
is  two-dimensional  with  fast-turn  around  times,  an  adaptation  algorithm  for  the  structured  grids 
used  for  the  flow  simulation  is  readily  available  and  easily  implemented,  and  that  the  dynamic 
mesh  terms  and  enforcement  of  the  Geometric  Conservation  Law  [Janus,  1989]  are  already  an 
integral  part  of  the  NS2D  software.  Although  different  flux  formulations  and  time  integration 
schemes  are  used  within  NS2D  and  Loci/CHEM,  for  the  cases  examined  here  they  evaluate  the 
same  equation  set  and  produce  comparable  results. 
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The  initiation  of  the  blast  event  is  characterized  by  interactions  between  a  planar  shock  wave, 
a  contact  discontinuity,  and  a  fan  of  rarefaction  waves,  as  shown  schematically  in  Figure  5.1. 
Although  similar  to  a  traditional  open-ended  shock  tube,  the  origin  in  Figure  5.1  represents  the 
center  line  of  the  high  pressure  region  which  acts  as  a  reflection  plane  for  the  shock  tube 
scenerio,  thus  resulting  in  a  blast  profile.  The  distance  Xo  is  the  half-width  of  the  high  pressure 
region  (as  shown  in  Figure  5.3).  In  order  to  assess  the  flow  solvers  ability  to  simulate  the  rapid 
release  of  energy  and  the  propagation  of  the  blast-wave,  ascertain  the  adequacy  of  the  grid 
resolution,  and  compare  the  NS2D  to  Loci/CHEM,  a  rigid  “flat”  wall,  simulation  of  the 
aforementioned  blast  was  undertaken.  To  this  end,  shown  in  Figure  5.2  is  a  comparison  of  the 
pressure  time  history  at  a  point  on  the  plate  from  the  NS2D  and  Loci/CHEM  simulations.  As 
seen,  excellent  agreement  is  observed  in  the  positive  phase,  negative  phase,  and  secondary  shock 
resolution.  The  slight  algorithmic  differences  between  NS2D  and  CHEM  can  explain  the 
difference  in  resolution  seen.  Further  analysis  is  ongoing  regarding  establishing  an  analytical 
standard  by  which  both  solutions  can  be  verified. 

The  initial  case  studied  was  the  interaction  of  a  blast  wave  with  a  circular-arc  “bump” 
geometry.  This  case  is  akin  to  that  utilized  by  Jaiman  et  al.  [Jaiman  et  al.,  2006]  to  assess  the 
accuracy  of  several  point-to-element  approaches  for  interdisciplinary  data  transfer.  For  a  rigid 
cylinder,  this  particular  case  has  been  numerically  studied  by  Ofengeim  and  Drikakis  [Ofengeim 
and  Drikakis  1997]  and  exhibits  significant  unsteady  flow  and  complicated  diffraction  of  shock 
waves  over  the  cylinder  during  the  blast.  Furthermore,  as  noted  in  the  literature  review,  for 
experimental  studies  with  a  single  or  a  small  number  of  particles,  it  has  been  difficult  to  arrive  at 
the  dependency  between  the  drag  coefficient  and  Reynolds  number.  Due  to  refraction,  reflection, 
and  diffraction,  the  geometry  and  amplitude  of  the  wave  fronts  are  altered  due  to  the  interaction 
of  the  shock  wave(s)  with  the  particle.  These  interaction  effects  become  significant  when  large 
debris/fragments  are  present. 

The  geometrical  description  and  initial  blast  conditions  are  shown  in  Figure  5.3;  only  a  half¬ 
plane  is  simulated  due  to  symmetry.  Approximate  blast  profiles  of  a  calculated  strength  and 
duration  can  be  created  by  adjusting  the  initial  high-pressure  and  temperature  within  the  blast 
zone  and  the  width  of  the  zone.  The  magnitude  of  the  blast  at  impact  on  the  cylinder  surface  can 
be  adjusted  via  the  blast  strength  and  standoff  distance  of  the  blast  initiation  zone.  The  conditions 
shown  in  this  report  are  for  a  half-width  (xo)  of  0.033Z,  centrally  located  upstream  a  distance  (xe) 
of  1.1  L,  where  L  is  the  distance  from  the  leading-edge  point  (B)  to  the  trailing-edge  point  (C). 
The  ambient  conditions  are  set  at  po=l atm  and  po=1.293£g/m3.  The  blast  zone  initial  pressure 
(stagnation  pressure)  is  set  such  that  p3/po=33.81  and  the  fluid  temperatures  are  uniform.  It 
should  be  noted  that  for  inviscid  flow  simulations,  the  conditions  for  imposing  symmetry  and  that 
representing  a  solid  body  are  equivalent  for  a  planar  surface.  Thus,  the  symmetry  line  could  also 
be  viewed  as  a  ground  plane.  The  structured  (quadrilateral)  grid  used  for  the  flow  simulation  in 
the  upper  domain  is  depicted  in  Figure  5.4. 

Modeling  of  the  solid  domain  is  accomplished  with  the  in-house  CSD  program  originally 
developed  by  Newman  [Newman  1997],  which  is  an  implicit  finite-element  code  that  performs 
nonlinear  (geometric  and  material  nonlinearity)  static  or  dynamic  analyses  and  includes 
numerous  element  types.  For  the  current  simulations,  however,  only  the  linear  elastic  plane 
stress  option  for  CSD  elements  is  utilized.  Furthermore,  within  this  software  the  Newmark-Beta 
algorithm  is  reformulated  for  the  time  integration  of  the  nonlinear  dynamic  response  such  that  it 
remains  consistent  with  the  flow  solver  during  FSI  simulations.  The  matched-mesh  structural 
discretization  of  the  elastic  solid  in  the  lower  domain  is  shown  in  Figure  5.4  along  with  the  CFD 
mesh.  The  coupled  simulations  assume  the  symmetry  line  of  the  fluid  domain  as  a  ground  plane 

103 


UNCLASSIFIED 


and  model  the  lower  half  as  an  elastic  body  under  plane  stress  conditions.  The  chosen  material 
properties  are  similar  to  those  used  for  solid  propellants  ( E=\QMPa ,  v=0.499,  p=\300kg/m3).  This 
particular  case,  albeit  having  limited  physical  relevance,  represents  an  excellent  scenario  by 
which  the  FSI  interdisciplinary  data  transfer  algorithms  may  be  evaluated.  Additionally,  it  may 
be  used  to  qualitatively  investigate  loosely  coupled  (two-way)  and  decoupled  (one-way)  response 
histories. 


Figure  5.1:  Characterization  of  blast  event  (x-t  diagram)  (Ofengeim  and  Drikakis,  1997). 


Figure  5.2:  Wall  pressuure  comparison  for  NS2D  vs  Loci/CHEM  (rigid  flat-plate). 
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Figure  5.3:  Geometry  used  for  blast  wave  interacting  with  a  cylinder  (Ofengeim  and 

Drikakis  1997). 


(b)  Close-up  of  the  leading  edge  of  the 
circular  arc  geometry. 


(a)  Near-field  view  of  the  disciplinary 
discretizations. 


Figure  5.4:  CFD  (red)  and  CSD  (green)  discretizations  used  for  blast- wave/bump  FSI. 
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Initially,  the  high  pressure  in  the  “blast  zone”  impulsively  loads  the  solid  media.  This  is 
followed  by  the  evolution  and  propagation  of  a  blast  wave  “pressure  profile”  which  impacts  the 
circular  arc  protrusion.  Figure  5.5  depicts  the  surface  pressure  distribution  at  selected  times  for 
the  (rigid)  CFD  simulations  which  are  used  for  the  one-way  FSI  coupling.  The  x-axis  represents 
the  horizontal  position  (measured  in  meters)  along  line  ABCD  in  Figure  5.3.  The  cylindrical 
surface  begins  at  0.0  m  (point  B)  and  ends  at  1 .0 m  (point  C).  Here,  the  pressures  from  the  CFD 
simulation  are  transferred  directly  to  the  deforming  CSD  surface  but  do  not  alter  the  blast-wave 
evolution.  Figure  5.6  depicts  color-shaded  pressure  fields  for  times  corresponding  to  when  the 
blast-wave  first  reaches  the  particular  points,  leading  edge  (LE),  top  (TOP),  and  trailing  edge 
(TE),  respectively.  Shown  in  Figure  5.7  are  the  decoupled  pressure  and  structural  response  for 
the  (LE),  (TOP),  and  (TE)  locations  of  the  circular  arc,  respectively.  For  this  one-way  coupling, 
significant  structural  response  only  occurs  once  the  blast  wave  reaches  each  location.  This  slow 
response  time  is  illustrated  clearly  in  Figure  5.8  as  the  blast- wave  reaches  the  leading  edge  of  the 
circular  arc.  A  comparison  of  surface  pressures  for  one-way  and  two-way  coupling,  as  shown  in 
Figure  5.9,  also  indicates  that  the  surface  deformations  have  little  influence  on  the  blast-wave  for 
the  conditions  chosen  for  this  demonstration. 

The  two-dimensional  test  case  presented  was  selected  because  it  sufficiently  represents  the 
physical  phenomena  found  in  a  blast  event.  Furthermore,  an  elastic  body  was  chosen  such  that 
the  structural  response  would  have  the  ability  to  interact  with  the  resulting  blast  waves 
propagating  over  and  impacting  an  obstacle.  The  decoupled  (one-way)  and  loosely  coupled  (two- 
way)  methodologies  where  both  simulated  and  compared  numerically.  Based  on  these 
simulations, ,  the  structural  response  for  the  selected  elastic  body  selected  was  slow  compared  to 
the  blast  event  and  only  local  deformations  were  present  that  did  not  significantly  alter  the  wave 
evolution.  These  preliminary  findings  are  consistent  with  observations  concerning  blast-wave 
(shock)-structure  interaction  found  in  the  literature.  This  would  not  be  the  case  for  thin  structures 
such  as  shells  which  would  have  faster  response  times  and  thus  have  large  non-local  deflections 
and  changes  in  surface  curvature.  The  development  of  the  2DFSI  will  enable  a  better 
understanding  of  limitations  of  the  one-way  coupling  approach. 


Figure  5.5:  Surface  pressure  distribution  at  selected  times  for  decoupled  (one-way) 
simulation  of  a  blast-wave/bump  FSI  interaction. 
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(a)  LE  impact,  t  ~  0.8ms. 


(b)  TOP  impact,  t  ~  1.48ms. 


(c)  TE  impact,  t  ~  2.65ms. 

Figure  5.6:  Shaded-surface  pressure  plots  of  blast-wave/bump  interaction. 
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(a)  Surface  pressure  and  deflections  at  the  LE. 


(b)  Surface  pressure  and  deflections  at  the  TOP. 
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(c)  Surface  pressure  and  deflections  at  the  TE. 

Figure  5.7:  Surface  pressure  and  deflections  for  decoupled  (one-way)  simulation  of  a  blast- 

wave/bump  interaction. 
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Figure  5.8:  Close-up  view  of  LE  surface  response  for  decoupled  (one-way)  simulation. 
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Figure  5.9:  Surface  pressures:  one-way  vs  two-way  coupling,  t  -  1.48ms. 
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5.2  Three-dimensional  blast-plate  FSI  simulation 

An  exemplary  experimental  setup,  by  which  preliminary  evaluation  and  validation  of  the 
developed  technology  may  be  made,  has  been  presented  by  Trzcinski  and  Cudzilo  [Trzcinski  and 
Cudzilo,  2004],  Their  study  presented  experimental  results  of  air  blast  waves  generated  by 
cylindrical  TNT  charges  of  various  weights  as  they  impacted  steel  and  laminated  fabric  plates. 
The  rectangular  plate,  having  a  side  length  of  0.5  m,  was  clamped  along  the  edges  as  shown  in 
Figure  5.10.  The  loaded  portion  of  the  plate  has  a  side  length  of  0.46m,  and  all  plates  have  a 
thickness  of  1  mm.  Experimental  plate  profile  final  shapes  were  recorded  for  three  TNT  charges 
(50g,  75g,  and  lOOg). 

For  the  three  dimensional  fluid- structure  interaction  simulations,  Loci/CHEM  was  used  with 
the  HLLE  inviscid  flux,  the  Barth  flux  limiter,  and  Runge-Kutta  second-order  time  integration. 
A  maximum  CFL  of  0.9  was  implemented  throughout  the  runs.  Although  the  charge  shape  was 
indicated  to  be  cylindrical  (as  in  a  “circular”  cylinder),  no  specific  dimensions  were  given  in  the 
report.  In  the  simulation,  the  charge  was  chosen  as  a  circular  cylinder  of  overall  square 
dimensions,  i.e.  cylinder  height  is  equal  to  the  cylinder  diameter.  For  computational  brevity,  the 
axisymmetry  (about  the  vertical  axis)  of  the  cylindrical  charge  was  exploited  in  addition  to  the 
quarter  plane  symmetry  of  the  entire  experimental  configuration.  In  figure  5.11,  the 
computational  geometry  and  mesh  cut-planes  for  this  configuration  are  shown  with  half-plane 
symmetry. 

The  solution  process  involved  performing  an  axisymmetric  simulation  of  the  charge 
detonation  on  a  uniform  mesh.  The  axisymmetric  mesh,  the  yellow  mesh  nearest  the  charge  in 
Figure  5.11,  has  an  order  of  magnitude  in  higher  resolution  (~0.5 mm  mesh  spacing)  than  the  fully 
3D  mesh  in  the  same  region.  The  charge  detonation  was  initiated  at  the  center  of  the  base  of  the 
charge,  and  the  axisymmetric  simulation  was  executed  for  ~0.05w5.  Limiting  the  simulation  to 
such  a  short  time  was  done  to  ensure  that  the  ambient  outer-boundary  condition  was  not 
disturbed.  The  solution  was  then  interpolated  to  the  computational  domain  of  the  fully  3D 
simulation.  The  solution  was  continued  for  an  additional  2 ms. 

The  50g  solutions  shown  in  Figure  5.12  are  of  the  logio  of  the  pressure  at  various  times  during 
the  blast  simulation.  Figure  5.12(a)  illustratessoon  after  the  solution  was  interpolated  from  the 
axisymmetric  detonation  simulation.  Figure  5.12(b)  shows  the  initial  impact  of  the  blast  wave 
with  the  plate  surface.  Although  the  view  is  from  above,  the  surface  pressures  on  the  bottom  of 
the  plate  are  depicted  in  the  figure.  Figure  5.12(c)  is  a  view  from  below  the  plate  illustrating  the 
reflected  blast  wave  (shown  in  the  colored-mesh  symmetry  plane).  Depicted  in  Figure  5.12(e)  is 
the  pressure  at  the  time  (-0.5 ms)  where  the  total  vertical  force  acting  on  the  plate  is  balanced. 
The  total  vertical  force  acting  on  the  plate  over  the  entire  time  history  of  the  simulations  for  each 
charge  weight  is  shown  in  Figure  5.13.  As  seen,  at  approximately  0.5 ms,  the  total  pressure  force 
begins  to  act  downwards  before  returning  to  atmospheric.  Since  the  blast  evolution  is  similar  for 
the  75g  and  lOOg  TNT,  pressure  plots  are  not  shown  for  brevity. 
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Figure  5.10:  Experimental  apparatus  for  blast  loading  on  plates  from  [Trzcinski  and 

Cudzilo  2004]. 


Figure  5.11:  Blast/plate  computational  geometry  (50g  TNT  cylindrical  charge). 
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(a)  t~0.06ms  (view  from  above  plate). 


(b)  t~0.15ms  (view  from  above  plate). 
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(c)  t~0.2ms  (view  from  under  plate). 


(d)  t~0.3ms  (view  from  above  plate). 
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(e)  t~0.5ms  (view  from  above  plate). 

Figure  5.12:  Computational  solution  for  50g  TNT. 


The  plate  surface  nodes  and  pressure  loads  from  the  Loci/CHEM  simulations  were  converted 
into  a  LS-DYNA  input  file  using  a  developed  utility  program.  This  program  reads  in  the 
extracted  data  at  the  user  defined  time  intervals  and  then  creates  the  appropriate  NODE, 
ELEMENTSHELL,  LOADSEGMENT  and  DEFINE  CURYE  cards  and  inserts  them  within 
the  input  file.  The  one-way  coupled  simulations  for  the  current  configuration  used  matched 
meshes  between  the  fluid  and  the  plate  structure.  However,  the  utility  can  easily  be  used  to 
convert  data  for  nested,  non-matched  meshes.  The  rectangular  plate  was  discretized  with  2500 
four-noded  Belytschko-Tsay  shell  elements  with  5-point  through  thickness  integration.  The  steel 
plate  is  modelled  using  the  plastic-kinematic  material  option  (MAT  003)  with  isotropic 
hardening  and  assumed  to  be  strain-rate  insensitive  (average  strain  rates  encountered  were  50  s'1). 
The  following  properties  for  steel  are  used:  density  p  =  7800  kg/m3,  elastic  modulus  E  =  196 
GPa,  tangent  modulus  Ep  =  0.2  GPa,  Poisson  ratio  v  =  0.3,  and  initial  yield  stress  oy  =  175  MPa. 
The  same  properties  were  utilized  by  Adamik  et  al.  [Adamik  et  al.,  2004]  in  a  LS-DYNA 
simulation  of  this  experiment. 

Shown  in  Figure  5.14  are  the  final  steel  plate  center  line  deflections  for  the  50g,  75g,  and  lOOg 
TNT  charges  from  the  experiment  of  Trzcinski  and  Cudzilo,  the  LS-DYNA  simulations  of 
Adamik  et  ah,  and  the  present  one-way  coupling  of  Loci/CHEM  and  LS-DYNA.  As  seen,  good 
agreement  with  the  experiment  is  observed.  The  lack  of  agreement  in  the  50g  TNT  case  in  Figure 
5.14a  is  most  likely  attributable  to  imperfections  in  the  experiment.  It  is  obvious  that  the  final 
deflections  in  the  experiment  were  not  symmetric,  whereas  for  the  75g  and  lOOg  TNT  cases  they 
were.  As  noted  in  Adamik  et  ah,  one  possible  imperfection  may  be  due  to  plate  clamping,  and  it 
was  found  during  simulations  that  varying  the  boundary  condition  treatment  at  the  clamped  edges 
had  a  significant  impact  on  the  plate  deflections.  Although  good  agreement  is  observed,  the  need 
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for  two-way  coupling  may  be  justified  as  the  charge  weight  is  increased.  This  is  evident  in  Figure 
5.14c  near  the  plate  center  where  drastic  increases  in  deflection  take  place. 

As  a  final  note,  the  results  from  the  experiment  are  comparable  to  the  LS-DYNA  simulation  of 
Adamik  et  ah;  however,  Adamik  et  al.  discussed  the  challenge  of  creating  an  initial  mesh  that  did 
not  require  mesh  rezoning  and  the  need  to  perform  the  simulation  in  distinct  phases.  The  need  for 
these  phases  were  primarily  driven  by  numerical  problems  associated  with  remote  regions 
(charge  and  contact  air  zone)  as  the  simulation  proceeded.  Adamik  et  al.  summarized  that  LS- 
DYNA  is  a  powerful  tool  for  simulating  such  blast  events,  but  they  conceded  that  considerable 
expertise  is  required  in  its  application.  On  the  other  hand,  the  one-way  coupling  between 
Loci/CHEM  and  LS-DYNA  was  easily  implemented  for  this  blast  wave  fluid- structure 
interaction. 


Figure  5.13:  Vertical  force  histories  for  the  50g,  75g,  and  lOOg  TNT  charges. 


115 


UNCLASSIFIED 


•  LS-DYNA  [Adamik  etal..  2004] 

*  Expe  rime  rrt  [T  rzcinski  and  Cudzib,  2004] 


•  LS-DYNA  [Adamk etal.,  2004] 

-  Experiment  [T rzcinski  and  Cudzilo,  2004] 


•  LS-DYNA  [Adamik  etal..  2004] 

-  Expe  rime  nt  [T  rzcinski  and  Cudzib,  2004] 


Figure  5.14:  Centerline  steel  plate  profile  for  different  charges. 


Although  no  experimental  data  existed,  the  current  blast-plate  was  utilized  to  ascertain  the 
influence  of  an  explosive  event  from  a  buried  charge.  To  this  end,  a  50g  TNT  charge  was  placed 
in  soil  (dry  sand)  at  the  same  standoff  distance  beneath  the  steel  plate  previously  discussed.  The 
soil  had  a  total  depth  of  0.15  m,  and  sat  upon  a  rigid  base  (reflection  plane).  The  depth  of  burial 
was  approximately  15  mm  beneath  the  top  soil  surface.  All  other  parameters  remained  the  same 
as  the  50g  TNT  air-blast  simulation.  The  total  vertical  force  acting  on  the  plate  for  both  the  50g 
TNT  air-blast  and  for  the  soil-blast  are  shown  in  Figure  5.15.  Here,  the  soil  confines  the  blast 
wave  and  thus  the  arrival  time  is  significantly  delayed.  Furthermore,  it  appears  in  Figure  5.15  that 
the  blast  wave  of  the  gas  arrives  at  a  lower  pressure,  followed  by  the  subsequent  impacting  of  the 
soil.  As  seen,  the  soil  imparts  a  much  greater  loading  on  the  plate.  The  combined  effect  of  the  gas 
and  soil  acts  on  the  plate  for  a  much  longer  duration.  However,  the  secondary  wave  (due  to 
reflections  within  the  initial  blasting  process)  appears  to  impact  the  plate  at  approximately  the 
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same  time  as  the  air-only  blast  (~1 .2  ms).  Figure  5.16a  depicts  the  initial  impact  of  the  blast 
wave,  soil  and  JWL  gas  on  the  plate,  illustrating  the  containing  and  focusing  effects  the  soil 
imparted  to  the  blast  process.  Figure  5.16b  depicts  the  final  stages  of  the  blast  interaction  with 
the  plate,  illustrating  the  dispersal  of  the  blast  products.  The  ultimate  result  of  simulating  this 
blast  in  soil  is  shown  in  Figure  5.17  for  the  final  steel  plate  center-line  deflections.  This  figure 
demonstrates  that  the  inclusion  of  soil  causes  the  blast  to  extend  over  a  much  larger  portion  of  the 
plate,  as  opposed  to  localized  central  deflections  seen  for  the  air-blast.  Additionally,  the  final 
central  point  deflection  is  actually  greater  than  those  of  the  75g  TNT  air-blast  (i.e.,  50%  more 
explosive). 


Figure  5.15:  Vertical  force  histories  for  the  50g  TNT  charges  in  air  and  soil. 
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(a)  Impact  with  steel  plate  (initial  stage) 


(b)  Impact  with  steel  plate  (final  stage) 


Figure  5.16:  Computational  solutions  for  the  50g  TNT  charge  in  soil  (dry  sand). 
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Figure  5.17:  Centerline  steel  plate  profile  for  50g  TNT  charges  in  air  and  soil. 
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Section  6:  Strategy  for  simulating  blast-vehicle 
interactions  for  complex  geometries 

Modeling  the  blast-soil  interaction  problem  is  a  particularly  challenging  problem  due  to  many 
factors.  Physically,  soil  can  be  considered  as  a  three  phase  mixture  composed  of  liquid  (water), 
gas  (air),  and  solids  (quartz,  organic  matter,  etc.).  When  considering  the  properties  of  a  soil  under 
a  high  strain-rate  event  such  as  explosive  blast,  the  soil  behavior  will  be  significantly  different 
than  in  low  strain-rate  events.  For  example,  in  typical  soil  tests,  the  liquids  and  gases  trapped  in 
the  pores  of  the  soil  will  be  free  to  migrate  and  so  will  become  part  of  the  measured  material 
strength  under  typical  experimental  compression  and  shear  tests.  Under  high  strain-rate  events, 
however,  the  liquid  and  gas  are  trapped  in  the  material  pores  and  do  not  have  sufficient  time  to 
migrate  through  the  material.  Thus,  for  explosive  events,  modeling  the  soil  as  a  mixture  equation- 
of-state  for  the  constituent  components  of  the  three  phase  material  provides  a  sound  starting  point 
for  developing  a  suitable  soil  model. 

In  the  interaction  between  the  soil  and  the  expanding  explosive  gases,  there  are  two  main 
regimes:  compression  and  tension.  In  the  initial  stages  of  the  explosive/soil  interaction,  the  soil  is 
essentially  behaving  as  a  fluid  undergoing  severe  deformation.  Our  present  simulations  for 
shallow  buried  mines  seem  to  confirm  this  with  very  good  crater  width  comparisons  for  early 
times  in  the  explosive  event.  As  time  proceeds  and  with  increasing  depth  of  burial,  we  find  that 
the  role  of  strength  begins  to  play  a  greater  role,  as  evidenced  by  over-predictions  of  crater 
widths  at  later  times  in  our  fluid  based  simulations.  However,  the  most  significant  challenge  in 
characterizing  the  behavior  of  soil  is  to  differentiate  its  response  to  tension  in  high  strain-rate 
events.  Under  tension,  soil  can  readily  loose  cohesion  and  become  an  aggregate  particulate  cloud 
of  soil  particles  and  soil  clumps.  The  effect  of  water  in  the  soil  must  also  be  accounted  for  in  the 
model.  Under  compression,  water  is  relatively  well-characterized  by  equations  of  state  such  as 
the  Tait  EoS  used  in  the  presently  reported  simulations.  Water  can  undergo  violent  phase  change 
(cavitation)  under  low  pressures.  These  low  pressure  conditions  are  more  likely  to  be  found  in 
high  moisture  content  soils  where  the  strength  of  the  relief  wave  caused  by  the  reflection  of  the 
blast  wave  off  of  the  soil-air  surface  can  produce  very  low  pressures.  Similarly,  on  the  outer  edge 
of  the  crater,  expansion  waves  can  produce  very  low  pressures.  For  most  soils  these  low 
pressures  will  likely  lead  to  water  cavitation  causing  soil  failure.  However,  unlike  simple  tearing 
of  the  soil,  cavitation  is  more  likely  to  produce  a  dense  cloud  of  small  particulates  that  are  highly 
energized  by  the  expanding  water  vapor.  The  present  models  that  were  used  in  this  study  and 
those  that  were  described  in  the  literature  do  not  model  the  effects  of  cavitation  of  the  water 
contained  within  the  soil. 

Understanding  the  mechanism  of  soil  failure  and  its  subsequent  transition  to  a  multiphase 
particulate  laden  flow  is  key  to  modeling  blast  soil  interactions  with  physical  reality.  There  are 
primarily  two  mechanisms  that  the  soil  plays  in  the  blast  loading:  1)  confinement  and  focusing 
effects,  and  2)  impulse  augmentation  due  to  solid  particle  impact  loading.  The  accurate  modeling 
of  both  of  these  effects  requires  an  accurate  modeling  of  the  transition  from  a  continuum 
treatment  of  the  soil  before  failure  to  a  multiphase  treatment  of  the  soil  after  failure.  We  note  that 
the  efficiency  of  momentum  transfer  between  the  soil  and  the  expanding  explosive  gases  is  most 
efficient  when  the  soil  fully  contains  the  explosive  gas.  As  the  soil  fails,  the  explosive  gases  will 
quickly  pass  around  the  soil  particles  to  the  surrounding  ambient  atmosphere.  Correctly  modeling 
this  process  will  have  a  major  effect  on  the  partition  of  blast  energy  that  is  deposited  to  the 
resulting  particulate  cloud  and  the  amount  of  blast  energy  that  remains  in  the  expanding  air  blast. 
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This  partitioning  of  energy  is  controlled  by  the  timing  of  the  failure  of  the  soil  and  the  subsequent 
expansion  of  the  debris  cloud.  Based  on  our  studies,  once  the  soil  becomes  free  particulates,  the 
coupling  between  the  explosive  gases  will  be  dominated  by  inviscid  and  volume  fraction  effects. 
In  regions  of  large  particle  volume  fraction  gradients,  significant  momentum  transfer  can  occur 
between  the  particles  and  the  expanding  gases.  In  addition,  when  relative  velocities  are  very 
large,  inviscid  effects  are  dominant.  However,  based  on  our  studies,  viscous  effects  can  be  safely 
neglected  for  particles  that  are  sand  grain-sized  or  larger.  As  a  result,  the  distribution  of  particle 
sizes  will  likely  play  only  a  small  role  in  the  accurate  prediction  of  blast  loads. 

Once  the  soil  has  failed  and  become  a  dispersed  phase  cloud  of  particles,  accurately  modeling 
the  transfer  of  momentum  from  the  expanding  cloud  of  detonation  products  will  be  essential  to 
accurately  modeling  the  vehicle  blast  loading.  In  part,  this  is  due  to  the  particles  being  more 
effective  at  transferring  momentum  to  the  solid  surface  than  to  a  gas.  In  the  momentum  transfer 
to  the  particles,  there  will  be  three  major  terms:  1)  viscous  drag,  2)  inviscid  wave  drag  (Mach 
number  corrections),  and  3)  pressure  action  due  to  gradients  in  volume  fractions.  The  momentum 
transfer  to  the  particulate  phase  will  either  be  coherent  (all  in  the  same  general  direction)  or 
randomly  distributed.  The  random  distribution  will  contribute  to  the  granular  temperature  of  the 
particulate  cloud,  while  the  coherent  transfer  will  contribute  to  the  bulk  motion.  Correctly 
modeling  the  balance  between  these  two  forms  of  momentum  transfer  will  play  an  important  role 
in  predicting  the  distribution  of  particle  trajectories  that  emanate  from  the  blast  crater. 

In  terms  of  the  physical  models,  a  key  concern  for  modeling  the  blast-soil  interaction 
accurately  will  be  correctly  modeling  the  process  of  soil  failure  and  subsequent  transfer  to  a 
dispersed  phase  model.  Since  the  solid  soil  material  will  be  very  effective  at  transferring  impulse 
to  the  vehicle,  the  timing  of  this  failure  will  be  critical  to  accurately  predicting  blast  loading:  A 
transition  that  is  too  early  will  under-predict  the  transfer  of  momentum  to  the  particulates,  while 
one  that  is  too  late  will  result  in  excessive  particle  velocities  and  a  subsequent  over-prediction  of 
blast  loading.  Immediately  after  failure,  the  large  volume  fractions  and  high  pressure  gradients 
will  yield  a  strong  but  abating  momentum  coupling  between  the  gas  and  particulate  phases. 
Currently,  the  main  challenge  with  creating  such  a  model  will  be  determining  the  appropriate 
phenomenological  models  for  soil  failure,  the  initial  energy  partition  between  soil  and  gas  upon 
failure,  the  energy  partition  between  coherent  kinetic  motion  of  the  particles  and  random  particle 
motion  represented  by  a  granular  temperature,  and  finally  the  model  for  the  momentum  coupling 
between  the  dispersed  and  gas  phases  after  soil  failure  is  complete. 

6.1  Numerical  modeling  of  soil-blast  interactions 

For  numerical  models  of  the  blast  soil  interaction  we  have  essentially  three  candidate 
methodologies:  1)  Lagrangian  methods,  2)  Eulerian  methods,  or  3)  meshless  methods. 
Lagrangian  methods  simulate  material  deformations  directly  by  tracking  the  deformation  of  the 
material  with  the  simulation  mesh.  These  methods  are  usually  based  on  a  finite-element 
formulation  and  are  robust  and  efficient  methods  for  computing  the  deformation  of  solids.  For 
highly  plastic  or  fluid  flow  problems,  they  are  subject  to  mesh  tangling  issues  that  can  cause 
accuracy  and  robustness  problems.  In  the  case  of  soil  in  close  proximity  to  an  explosive  event, 
the  soil  undergoes,  for  a  period  of  time,  severe  deformations  that  can  present  challenges  for  a 
Lagrangian  method.  In  addition,  the  simulation  of  the  soil  using  a  Lagrangian  method  will  require 
some  sort  of  fluid-structure-interaction  to  account  for  the  interaction  of  the  soil  with  the 
expanding  explosive  gases. 
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An  alternative  to  Lagrangian  methods  are  Eulerian  methods  that  keep  the  mesh  fixed  and 
allow  material  to  move  through  the  mesh.  These  methods  can  be  implemented  using  either  finite- 
element  or  finite-volume  formulations  and  solve  the  mesh  tangling  issue.  In  addition,  it  is 
possible  to  formulate  an  Eulerian  solver  such  that  the  soil  and  explosive  are  treated  in  one  unified 
model  whereby  issues  such  as  conservation  of  mass,  energy,  and  momentum  are  trivially 
satisfied.  One  downside  of  the  Eulerian  techniques  is  that  they  have  less  fidelity  in  capturing 
material  interface  dynamics.  However,  for  blast-soil  interactions,  the  soil  interface  is  unlikely  to 
remain  sharply  defined,  and  so  such  a  compromise  is  reasonable. 

A  third  method  that  can  be  highly  robust  for  simulating  dynamically  evolving  interfaces  are 
meshless  methods  such  as  Smooth  Particle  Hydrodynamics  (SPH).  In  these  methods,  the  material 
is  represented  by  a  collection  of  "particles"  that  can  evolve  in  time.  Particles  interact  through  the 
use  of  interpolation  functions  that  average  particles  within  a  given  smoothing  length  (which  can 
be  adjusted  based  on  local  particle  densities).  It  should  be  noted  that  such  a  scheme  is  still  a 
continuum  simulation  whereby  the  "particles"  do  not  actually  represent  particulates  as  found  in  a 
multiphase  flow,  but  rather  are  artifacts  of  this  numerical  discretization  methodology:  the 
underlying  numerical  model  is  still  solving  equations  of  continuum  mechanics.  These  SPH 
methods  are  more  expensive  than  either  Lagrangian  or  Eulerian  methods  and  may  provide  a  more 
natural  way  to  transition  from  a  continuum  to  multiphase  representation.  However,  it  would  seem 
that  the  cost  of  these  methods  are  not  justified,  as  we  can  always  combine  the  Lagrangian  or 
Eulerian  methods  with  multiphase  models  which  do  not  require  the  expensive  smoothing  function 
and  thus  would  be  much  more  efficient.  We  would  note  that  recent  work  in  combining 
Lagrangian  and  SPH  methods  [Johnson  et  al.,  2008]  show  the  potential  richness  of  using  SPH  to 
couple  soil  and  blast  simulations.  We  do  note  that  while  this  paper  does  use  SPH  particles  to 
improve  robustness  of  Lagrangian  methods,  they  do  not  consider  the  soil  transition  to  a 
multiphase  flow  problem. 

The  modeling  capabilities  that  have  been  developed  thus  far  include  an  Eulerian  finite-volume 
solver  that  is  fully  conservative  for  the  simulation  of  blast-soil  interactions.  This  model  currently 
does  not  include  material  strength  models,  nor  does  it  consider  the  failure  of  the  soil  and  the 
subsequent  transition  to  accurate  multiphase  models.  In  order  to  improve  the  modeling  of  blast 
loads  due  buried  explosive  devices,  we  suggest  that  our  current  model  should  be  extended  to 
include  material  strength.  In  addition  we  need  to  extend  our  present  model  such  that  we  couple 
soil  failure  to  a  multiphase  model  of  the  resulting  debris  cloud. 

Once  the  basic  model  is  in  place  that  can  accurately  capture  the  blast-soil  interaction  and  the 
multiphase  debris  cloud,  the  next  important  focus  is  to  determine  the  appropriate  model  for  soil 
failure:  What  are  the  mechanisms  for  soil  failure?  What  role  does  cavitation  play  in  soil  failure? 
How  should  cavitation  be  modeled?  When  failure  happens,  how  important  is  capturing  the 
transfer  of  energy  to  the  granular  kinetic  energy  of  the  particulate  cloud?  How  is  this  partition  of 
blast  energies  accurately  modeled? 

6.2  Numerical  modeling  of  blast-structural  interactions 

The  accurate,  robust,  and  efficient  simulation  of  fluid-structural  interactions  for  complex 
geometries  involving  even  modest  physical  models  represents  a  significant  challenge  as  well.  The 
phenomena  associated  with  explosive  events,  and  the  resulting  interactions  between  the  various 
materials  and  states,  along  with  the  potential  for  fragmentation  and  fracture,  make  this  task  even 
more  daunting.  This  problem  becomes  more  tractable  when  high-fidelity  simulations  codes  are 
available  to  model  the  relevant  physics  associated  with  the  various  disciplinary  domains,  and  then 
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compatibility/equilibrium  is  enforced  at  the  interfacial  boundaries.  However,  enforcement  of  this 
compatibility  represents  yet  another  challenge  towards  accuracy,  robustness  and  efficiency.  To 
this  end,  numerous  methods  have  been  developed  and  devised  to  tackle  these  issues,  with  each 
having  specific  advantages  and  disadvantages.  Although  a  comprehensive  review  of  these 
methods  will  not  be  presented  within  this  report,  some  of  the  various  attributes  shall  be  discussed. 
Furthermore,  many  common  features  exist  between  modeling  the  soil-blast  and  the  blast- 
structural  interactions,  and  can  possibly  be  exploited.  These  discussions  shall  put  into  context  the 
candidate  methodologies  for  simulating  blast-vehicle  interactions  using  Loci/CHEM  and  LS- 
DYNA. 

As  presented  in  Section  6.1  for  soil-blast  interactions,  it  is  well  known  that  when  using  a 
purely  Lagrangian  approach  that  mesh  tangling  becomes  a  significant  issue  when  large 
deformations  are  encountered  and  create  the  need  for  costly  remeshing,  potential  introduction  of 
errors,  and  lack  of  robustness  when  dealing  with  topological  changes.  In  Eulerian  approaches 
numerical  diffusion  (due  to  advection  of  the  solid  stress  or  strain)  may  result  in  substantial 
inaccuracies  and  additionally  has  the  burden  of  locating  material  interfaces.  Many  techniques 
have  been  developed  in  an  attempt  to  mitigate  these  deficiencies  for  specific  problem  classes. 
One  such  approach,  for  example,  is  the  material  point  method  (MPM)  extended  for  multiphase 
flows  [Zhang  et  al.,  2008],  That  technique  was  developed  to  use  both  an  Eulerian  mesh  and 
Lagrangian  material  points.  The  velocity  of  each  phase  is  computed  on  the  Eulerian  mesh,  thus 
eliminating  mesh  tangling  problems.  The  deformation  (and  thus  stress  history)  of  the  material 
points  are  tracked  by  the  Lagrangian  nodes.  That  particular  method  avoids  most  of  the  pitfalls 
associated  with  the  Eulerian  and  Lagrangian  formulations,  but  retains  their  respective  advantages. 
However,  the  computational  costs  are  higher  and  there  is  need  to  interpolate  information  between 
the  mesh  nodes  and  material  points  at  each  time  step.  Additionally,  this  method  suffers  from 
errors  (particularly  at  sharp  interfaces)  due  to  the  inability  to  satisfy  the  continuity  equation  in 
terms  of  the  volume  fractions  identically,  and  is  only  satisfied  in  the  sense  of  a  weak  solution. 

At  this  point,  the  issue  related  to  transition  within  the  soil-blast  interaction  and  the  blast- 
structural  interactions  becomes  more  apparent.  In  the  early  stages  of  the  soil-blast  evolution,  a 
continuous  multiphase  flow  is  appropriate  to  model  the  large  fragments  associated  with  the 
interaction.  As  the  soil  breaks  apart  into  smaller  fragments,  it  transitions  to  a  dispersed 
multiphase  flow.  For  the  fluid-vehicle  structure  interaction,  the  simulation  is  similar  to  that  of  the 
early  stages  of  the  soil-blast  model,  and  may  be  represented  with  the  same  continuous  multiphase 
flow  approximation. 

As  is  apparent  from  the  discussions  thus  far,  we  wish  to  develop  a  method  that  avoids  the 
pitfalls  of  the  various  approaches,  may  be  readily  implemented  within  Loci/CHEM,  and  is 
ultimately  two-way  coupled  with  LS-DYNA.  Although  many  details  require  further  research,  one 
candidate  methodology  is  to  use  Loci/CHEM  as  an  Eulerian  predictor  followed  by  using  LS- 
DYNA  as  a  Lagrangian  corrector  step.  For  such  an  algorithm,  implementation  of  a  strength 
model  into  Loci/CHEM  would  be  required.  However,  development  of  strength  models  in  the 
context  of  a  Godunov  type  solver  similar  to  the  algorithm  employed  within  Loci/CHEM  has  been 
developed  [Miller  and  Colella,  2001],  Miller  and  Colella  reformulate  the  equations  of  solid 
mechanics  as  a  first-order  system  of  hyperbolic  equations  and  solve  the  system  via  an  explicit 
second-order  Godunov  method.  In  particular  that  solution  algorithm  is  a  predictor-corrector 
scheme  that  solves  the  conservative  fluxes  representing  elastic  flow  as  the  predictor  followed  by 
a  corrector  step  accounting  for  plasticity  within  a  source  term.  This  method  conserves  mass, 
momentum,  and  energy  exactly,  but  does  not  conserve  the  deformation  gradient.  Within  the 
Eulerian  grid  approaches,  the  immersed  boundary  is  usually  established  and  tracked  via  some 
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form  of  interface  markers  or  advected  as  a  solution  variable  such  as  a  level-set.  Utilizing  LS- 
DYNA  as  the  Lagrangian  corrector  step,  the  update  represents  the  true  location  of  the  interface 
within  the  Eulerian  grid,  and  conservation  of  the  deformation  gradient  is  no  longer  an  issue.  The 
Eulerian  solution  would  then  require  subsequent  correcting  to  establish  consistency  between  the 
two  formulations.  At  the  interface  this  requires  that  kinematic  and  force  equilibrium  is 
maintained. 

To  achieve  the  aforementioned  coupling,  several  challenges  still  require  investigation.  For 
example,  two-way  coupling  with  LS-DYNA  appears  achievable  through  the  external  user  defined 
loading  routines,  especially  for  an  ALE  approach.  However,  within  the  corrector  step,  the  transfer 
of  loading  information  may  require  additional  solution  control  options  such  that  equilibrium  is 
achieved  during  each  iteration.  Since  both  solution  algorithms  are  explicit,  and  thus  limiting  the 
propagation  of  information  within  each  domain,  local  strategies  should  be  investigated  to  enforce 
equilibrium.  In  the  presence  of  fracture,  any  void  or  crack  within  the  structure  will  not  be 
resolved  at  sufficient  resolution  to  accurately  model  the  resulting  fluid  behavior.  Detailed  studies, 
and  potential  use  of  phenomenological  models,  need  to  be  explored  in  such  instances. 
Additionally  the  use  of  local  mesh  adaptation  should  be  considered  when  phenomenological 
models  are  inappropriate. 

Finally,  we  would  note  that  for  blasts  that  are  a  moderate  distance  from  the  explosive,  it  is 
reasonable  for  the  gas  phase  coupling  to  be  performed  using  a  one-way  coupling  approach.  The 
validity  of  using  a  one-way  coupling  approach  is  largely  due  to  the  dramatic  differences  between 
the  gas  and  structure  material  densities:  Even  with  perfect  momentum  transfer,  a  factor  of  100 
difference  in  densities  would  imply  a  similar  factor  difference  in  the  resulting  material  velocities. 
For  the  gas  phase  coupling,  it  is  possible  to  estimate  the  error  involved  in  making  a  one-way 
coupling  assumption.  Such  a  tool  would  allow  an  engineer  to  first  perform  a  one-way  coupled 
simulation  and  then  assess  from  this  simulation  if  a  more  expensive  two-way  coupling  approach 
is  needed.  We  would  also  note  that  while  the  gas  phase  densities  drop  dramatically  with  distance 
from  the  explosive  event,  the  multiphase  problem  does  not.  Therefore,  a  compromise  1 1/2  way 
coupling  may  be  reasonable  in  many  cases.  In  this  mode,  the  gas  phase  coupling  would  be 
treated  one-way,  while  the  particulate  loading  would  be  aware  of  the  dynamic  changes  of  the  LS- 
DYNA  simulation  of  the  vehicle.  If  these  techniques  were  combined  with  a  robust  tool  for 
estimating  approximation  errors,  engineers  could  quickly  turn  around  risk  assesments  and 
determine  if  the  approximations  yielded  too  large  an  increase  in  uncertainties  associated  with 
predictions. 
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